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life  and/or  property.  Extensive  wind-tunnel  testing  and  radio-controlled 
flight  testing  has  been  done  over  the  last  twenty  years  to  gain  a  better 
understanding  of  the  dynamic  instabilities  at  high  angles-of-attack.  A 
basic  problem  has  existed  in  interpreting  these  data  and  in  making  predic¬ 
tions  of  aircraft  dynamic  behavior  so  as  to  achieve  close  agreement  with 
flight  test  data. 

Aircraft  dynamic  behavior  at  high  angles-of-attack  is  highly  nonlinear 
and  in  the  past  there  has  been  a  lack  of  suitable  techniques  for  analyzing 
the  global  behavior  of  nonlinear  systems.  Under  a  previous  project  with 
the  Office  of  Naval  Research,  Scientific  Systems,  Inc.  has  developed  a  new 
approach  based  on  Bifurcation  Analysis  and  £atastrophe  Theory  Methodology 
(BACTM) .  The  approach  has  been  applied  to  specific  jump  and  iTmit  cycle 
behavior  such  as  roll -coupling,  pitch  up,  post-stall  departure,  divergence, 
spin  entry,  developed  erect  spin,  and  spin  prevention  and  recovery.  The 
aircraft  used  for  the  study  of  spin  motions  was  selected  because  of  the 
completeness  of  the  aero  data  in  the  spin  flight  regimes,  and  because  it 
is  representative  of  modern  fighters.  This  model  was  also  used  for  studies 
of  non-spin,  high  angle-of-attack  behavior. 

Under  this  project,  the  full  six  DOF  aircraft  model  was  implemented, 
and  used  not  only  for  the  above  studies,  but  also  for  several  new  de¬ 
velopments  in  the  BACTM  methodology.  The  new  developments  are  basically 
in  the  area  of  generalizing  and  improving  the  numerical  techniques  for 
computing  equilibrium  and  bifurcation  surfaces,  in  expanding  the  com¬ 
prehensiveness  of  the  physical  model  and  environment  and  in  the  study  of 
chaotic  motions.  ' 

The  work  on  this  project  has  centered  around  the  application  of 
BACTM  to  study  the  spin  characteristics  of  a  "variable  sweep"  fighter 
aircraft.  The  aerodynamic  data  for  this  model  roughly  corresponds  to 
experimental  data  for  the  F-1I1,  although  modifications  in  some  of  the 
numbers,  particularly  Cn,  are  required  to  make  simulation  results  agree 

with  flight  test  data.  We  have  designated  this  simulation  model  as 
Aircraft  F. 

Spin  behavior  is  typically  a  post-stall  phenomenon,  and  is  character¬ 
ized  by  angles-of-attack  much  in  excess  of  the  stall  value  of  angle-of- 
attack.  It  Is  also  possible  that  spin  conditions  will  follow  a  noil 
departure  motion.  A  certain  type  of  spin,  the  erect  flat  spin,  has 
been  given  particular  emphasis  in  this  work  effort.  This  spin  is  featured 
by  values  of  a  (angle-of-attack)  in  the  75-85  degree  ranges;  a  vertical 
body  rotation  rate,  which  is  also  constant  over  time,  and  center  of  mass 
motion  which  is  basically  helical,  with  the  axis  parallel  to  local 
vertical;  and  a  noticeably  prominent  yaw  rate. 
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CHAPTER  1 


Introduction  and  Summary 

1.1  Scope  of  Work  Effort 

High  angle-of-attack  phenomena  have  been  of  interest  to  aerodynamicists, 
aircraft  designers,  pilots  and  control  system  analysts  ever  since  the  advent 
of  modern  high  performance  aircraft.  Due  to  the  concentration  of  inertia 
along  the  fuselage,  the  modern  jet  fighters  are  highly  susceptible  to 
post-stall  departures  and  spin.  In  spite  of  extensive  design  effort,  modern 
aircraft  still  inadvertently  enter  spins  which  sometimes  result  in  loss 
of  life  and/or  property.  Extensive  wind-tunnel  testing  and  radio-controlled 
flight  testing  has  been  done  over  the  last  twenty  years  to  gain  better 
understanding  of  the  dynamic  instabilities  at  high  angles-of-attack.  A 
basic  problem  has  existed  in  interpreting  these  data  and  in  making  predic¬ 
tions  of  aircraft  dynamic  behavior  so  as  to  achieve  close  agreement  with 
fl ight  test  data. 

Aircraft  dynamic  behavior  at  high  angles-of-attack  is  highly  nonlinear 
and  in  the  past  there  has  been  a  lack  of  suitable  techniques  for  analyzing 
the  global  behavior  of  nonlinear  systems.  Under  a  previous  project  with 
the  Office  of  Naval  Research,  Scientific  Systems,  Inc.  has  developed  a  new 
approach  based  on  Bifurcation  Analysis  and  Catastrophe  Theory  Methodology 
(BACTM).  The  approach  has  been  applied  to  specific  jump  and  limit  cycle 
behavior  such  as  roll-coupling,  pitch  up,  post-stall  departure,  divergence, 
apin  entry,  developed  erect  spin,  and  spin  prevention  and  recovery.  The 
aircraft  used  for  the  study  of  spin  motions  was  selected  because  of  the 


2 


completeness  of  the  aero  data  in  the  spin  flight  regimes,  and  because  it 
is  representative  of  modern  fighters.  This  model  was  also  used  fir  studies 
of  non-spin,  high  angle-of-attack  behavior. 

Under  this  project,  the  full  six  DOF  aircraft  model  was  implemented, 
and  used  not  only  for  the  above  studies,  but  also  for  several  new  de¬ 
velopments  in  the  BACTM  methodology.  The  new  developments  are  basically 
in  the  area  of  generalizing  and  improving  the  numerical  techniques  for 
computing  equilibrium  and  bifurcation  surfaces,  in  expanding  the  com¬ 
prehensiveness  of  the  physical  model  and  environment  and  in  the  study  of 
chaotic  motions. 

The  work  on  this  project  has  centered  around  the  application  of 
BACTM  to  study  the  spin  characteristics  of  a  "variable  sweep"  fighter 
aircraft.  The  aerodynamic  data  for  this  model  roughly  corresponds  to 
experimental  data  for  the  F-lll,  although  modifications  in  some  of  the 
numbers,  particularly  C^,  are  required  to  make  simulation  results  agree 
with  flight  test  data.  We  have  designated  this  simulation  model  as 
Aircraft  F. 

Spin  behavior  is  typically  a  post-stall  phenomenon,  and  is  character¬ 
ized  by  angles-of-attack  much  in  excess  of  the  stall  value  of  angle-of- 
attack.  It  is  also  possible  that  spin  conditions  will  follow  a  roll 
departure  motion.  A  certain  type  of  spin,  the  erect  flat  spin,  has 
been  given  particular  emphasis  in  this  work  effort.  This  spin  is  featured 
by  values  of  a  (angle-of-attack)  in  the  75-85  degree  ranges;  a  vertical 
body  rotation  rate,  which  is  also  constant  over  time,  and  center  of  mass 
motion  which  is  basically  helical,  with  the  axis  parallel  to  local 
vertical;  and  a  noticeably  prominent  yaw  rate. 
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The  study  of  spin  behavior  begins  with  an  analysis  of  the  types  and 
nature  of  equilibrium  spin  conditions.  Because  gravity  plays  a  role, 
the  basic  system  of  equations  is  eighth  order  (the  six  force-moment  equa¬ 
tions,  plus  two  kinematical  equations  for  pitch  and  roll  angles).  Gravity 
is  not  a  significant  factor  in  the  so-called  "roll-coupling"  flight  regime, 
studied  previously  (Mehra  et  al .  (1977)).  In  the  case  of  spin  equilibrium 
conditions,  the  presence  of  a  non-zero  gravity  term  causes  the  roll  and 
pitch  angles  to  enter  the  basic  sixth-order  system. 

The  highly  nonlinear  nature  of  flat  spin,  and  the  extreme  values  of 
state  variables  which  typify  it,  require  that  the  aerodynamic  data  extend 
over  values  of  a  and  sideslip  (6)  which  are  well  beyond  the  ranges  of 
readily  accessible  data.  The  data  we  have  used  here  for  aircraft  F  were 
available  in  tabular  form,  and  do  encompass  the  necessary  ranges  (Moore 
et  al.  (1971)).  Spline  function  polynomials  were  used  to  model  this  aero 
data  because  these  functions  are  continuous  at  all  interior  points  in¬ 
cluding  certain  derivatives.*  Spline  functions  can  therefore  give  accurate 
results  over  all  points  in  the  region  with  the  accuracy  needed  to  insure 
efficient  numerical  solution  of  the  equilibrium  and  bifurcation  surfaces. 

There  is  also  flexibility,  in  that  numerical  techniques  which  utilize 
analytical  expressions  for  the  derivatives  of  the  aero  coefficients,  can 
effectively  utilize  the  spline  approximation.  Our  results  have  confirmed 
the  soundness  of  this  choice. 

A  final  note  on  the  simulation  model.  The  controls  chosen  were  the 
standard  aerosurface  controls,  aileron,  elevator  and  rudder  deflections. 

Thrust  is  not  used  as  a  control  explicitly.  This  is  not  an  undue  restriction. 


♦This  is  true  up  to  second  order  derivatives  when  cubic  or  bi-cubic  splines 
are  used. 
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since  thrust  is  generally  maximum  during  high-a  maneuvering  or  is  reduced 
to  idle  during  spin. 

Equilibrium  surfaces  for  spin  conditions  were  computed  with  satis¬ 
factory  accuracy  and  efficiency  by  means  of  a  parametric  continuation  pro¬ 
cedure  based  on  the  methods  of  Davidenko  (1953)  and  Lahaye  (1934).  In 
its  basic  form,  this  procedure  solves  a  system  of  nonlinear  algebraic 
equations  by  varying  a  parameter  from  a  value  for  which  the  unknowns  are 
readily  determined  to  the  actual  value  of  the  basic  system.  In  our  appli¬ 
cation,  the  parameter  is  one  of  the  aerosurface  controls  and  the  unknowns 
are  the  eight  state  variables.  The  starting  point  is  determined  by  a 
Newton-Raphson  scheme,  with  initial  guesses  for  the  state  and  control  at 
values  which  correspond  to  expected  spin  situations.  In  aircraft  F,  for 
example,  the  equilibrium  pitch  and  roll  rates  are  about  20  deg/sec,  and 
yaw  rate  is  roughly  10  times  as  large.  Velocity  is  about  450  feet  per 
sec,  angle-of-attack  about  83  deg.  Sideslip,  pitch  and  roll  equilibrium 
angles  are  5  deg  or  less.  It  has  been  verified  that  the  equilibrium  angular 
velocity  is  vertical.  The  continuation  parameter,  say  rudder  deflection, 
is  then  extended  over  its  range  from  this  starting  point. 

The  particular  continuation  technique  employed  here  is  principally  an 
amalgamation  of  methods  proposed  by  Klopfenstein  (1961),  Keller  (1977), 
and  Kubicek  (1976).  These  methods  arose  out  of  the  necessity  of  dealing  with 
various  singularities  which  typify  nonlinear  equilibrium  surfaces.  The  most 
important  of  these  singularities  are  limit  points  and  bifurcation  points. 
Kubicek,  Keller  and  Klopfenstein  have  added  an  arclength  parameter,  making 
it  the  independent  variable,  to  eliminate  the  limit  point  singularity.  Kubi¬ 
cek  and  Klopfenstein  use  a  purely  Euclidean  arclength  parameter,  while  Keller 
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introduces  a  family  of  "pseudo"  arclength  parameters  which  help  to  allow  one  to 
solve  for  the  slopes  of  all  of  the  solution  curves  which  pass  through  a  bifur¬ 
cation  point.  In  this  way,  the  bifurcation  points  may  be  isolated. 

Continuation  proceeds  on  the  original  branch  using  the  knowledge  of  its 
slope  and  "jumping"  over  the  actual  bifurcation  point.  Then,  one  returns 
to  the  bifurcation  point  and  begins  continuation  along  the  secondary  branches. 

In  addition,  Keller's  scheme  may  be  extended  to  function  spaces,  so  that 
differential  systems  such  as  Two-Point  Boundary  Value  Problems  may  also  be 
solved.  The  methods  of  Kubicek  and  Keller  enable  the  computation  of  the 
complete  equilibrium  or  bifurcation  surface  with  just  one  computer  run, 
and  are  typically  more  "robust"  than  Klopfenstein 1 s  method. 

The  eigenvalue  analysis  of  the  equilibrium  surfaces  for  aircraft  F 
show  regions  of  jump,  limit  cycle,  hysteresis  and  other  phenomena  similar 
to  those  found  for  aircrafts  A,  B  and  H  investigated  earlier  (Mehra  et  al . 
(1977)).  In  the  case  of  the  spin  phenomena,  the  magnitudes  of  the  jumps 
are  typically  smaller,  though  there  is  indication  that  these  jumps  would 
go  from  flat  spin  to  an  intermediate  spin  (a  ^70°)  to  steep  spin  (a -45°). 
Additionally,  jumps  to  limit  cycles  or  oscillatory  spins  are  also  present. 

The  aircraft  F  model  has  also  been  used  to  generate  equilibrium  sur¬ 
faces  in  non-spin  regions  prior  to  departure.  A  major  consideration 
of  such  non-spin  regimes  is  that  roll  and  pitch  angles  generally  do  not 
have  equilibrium  (steady  state)  values,  so  that  these  variables  must  be 
decoupled  from  the  basic  system.  This  is  done  by  neglecting  gravity 
effects.  The  results  are  similar  to  those  obtained  with  Drevious  models 
(e.g.,  aircraft  H).  However,  since  angle-of-attack  data  were  available 
only  to  -10  degrees,  roll  daparture  studies  with  aircraft  F  were  somewhat 
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limited.  A  study  of  the  nature  of  equilibrium  surfaces  in  the  "transition" 
region  between  roll  departure  and  flat  spin  has  also  been  made. 

In  what  may  prove  to  be  a  very  useful  result  in  the  analysis  of 
spin  motions,  it  was  found  that  equilibrium  surfaces  in  the  spin  regime, 
which  have  all  of  the  features  of  the  standard  spin  equilibrium  surfaces, 
can  be  generated  using  the  lower-dimensional  "non-spin"  equilibrium 
system.  In  this  system,  gravity  is  assumed  to  be  zero.  It  should  be 
emphasized  that  the  numerical  results  are  often  different,  but  the  shape 
of  the  curves  is  quite  similar.  This  approximation  has  been  made  only  in 
the  study  of  flat,  developed  spins.  These  spins  are  characterized  by 
high  «Spj^  values,  and  low  spin  equilibrium  pitch  angles  (0SpjN)’  With 
this  situation,  gravity  terms  in  the  dynamic  equation  for  a  are  small. 
Finally,  for  many  of  the  spin  conditions,  changing  the  value  of  V, 
the  velocity  magnitude,  by  30%  had  a  greater  effect  on  the  equilibrium 
curves  than  did  eliminating  gravity. 

Time  history  runs  of  spin  conditions  for  aircraft  F  have  been  made, 
and  results  confirm  the  predictions  of  the  equilibrium  surfaces.  Those 
runs  which  begin  in  the  developed  flat  spin  condition  follow  quite  accurately 
the  results  predicted  by  the  spin  equilibrium  curves.  Runs  which  attempt 
to  achieve  flat  spin  from  non-spin  conditions  were  also  made.  It  has 
been  found  that,  as  reported  in  Bihrle  (1976),  ensuing  motion  in  spin 
regions  is  highly  sensitive  to  both  the  initial  conditions  and  the  control 
sequencing.  Also,  we  discovered  that  it  is  much  easier,  for  the  given  simu¬ 
lation  model,  to  achieve  equilibrium  spin  from  "trim"  flight  conditions 
(controls  neutral)  when  the  velocity  magnitude  is  held  fixed.  This  is 


equivalent  to  both  adjusting  thrust  magnitude  and  vectoring  the  thrust, 
to  keep  it  aligned  with  the  current  velocity  vector  and  keeping  con¬ 
stant  magnitude.  This  is  an  approximation,  but  it  does  obviate  the 
need  at  the  current  time  to  become  concerned  with  the  role  of  vehicle 
thrust  in  post  stall  and  spin  entry  conditions.  Such  a  concern  can 

be  more  readily  dealt  with  when  more  comprehensive  models  such  as  the  F-4 
are  implemented.  Similarly,  spin  recovery  simulations  have  been  made, 
and  the  results  again  indicate  that  the  aircraft  F  model  is  highly  sensitive 
to  the  recovery  control  sequence.  We  expect  to  develop  a  systematic  re¬ 
covery  methodology  once  the  spin  bifurcation  surfaces  are  completed. 
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1.2  Previous  Work 

Most  of  the  theoretical  and  experimental  results  relating  to  the 
characteristics  of  spin  motion  has  been  performed  at  the  NASA/Langley 
Research  Center.  The  majority  of  this  work  has  been  experimental  in 
nature.  Theoretical,  or  analytical,  results  have  been  hampered  by  two 
factors:  1)  the  highly  nonlinear  nature  of  the  spin  regime,  and  2)  the 
difficulty  both  of  obtaining  wind  tunnel  data  which  are  relevant  to  spin 
motions  and  of  effectively  correlating  these  data  with  the  actual  air¬ 
craft's  performance. 

Klinar  and  Grantham  (1959)  used  traditional  linearized  analysis  tech¬ 
niques  to  study  flat,  steady  spin  behavior.  However,  similar  efforts  done 
previously  have  been  limited  to  reliance  on  the  limited  conventional 
static  and  forced  oscillation  aerodynamic  data.  These  data  do  not  always 
represent  adequately  the  highly  complex  flow  phenomena  associated  with 
flight  in  these  stall/spin  regions.  Consequently,  the  wind-tunnel  tech¬ 
niques  were  expanded.  Neihouse  et  al.  (1957,  1960)  report  on  the  develop¬ 
ment  of  a  rotary  balance  mechanism  by  which  a  model  is  spun  freely  about 
selected  spin  axes,  over  wide  ranges  of  angle-of-attack.  They  also  dis¬ 
covered  that  differences  between  the  model  results  and  those  of  the  air¬ 
craft  became  more  pronounced  with  current  high-speed  designs.  The  dif¬ 
ferences  were  felt  to  be  due  to  such  factors  as  possible  aerodynamic 
scale  effects  (or  Reynolds  number  effects)  and  variations  in  testing 
techniques  between  airplanes  and  free-spinning- tunnel  models. 

Analysis  made  by  Scher  and  Anglin  (1959)  further  determined  that 
different  kinds  of  spins  are  entered  depending  upon  whether  the  aircraft 
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was  initially  in  a  trinmed,  level  flight  condition  or  spin  entry  was 
achieved  from  an  applied  yaw  rotation  which  simulates  the  spin-tunnel 
mode.  Using  BACTM  on  our  aircraft  F  model  (Chapter  III),  we  obtain 
similar  results.  Their  results  again  emphasize  the  importance  of  Reynolds 
number  effects  as  well  as  tunnel  test  techniques.  More  recently,  Bihrle 
(1974,  1976)  corroborates  the  large  effect  of  Reynolds  number  on  spin 
aerodynamic  characteristics.  In  addition,  he  recognizes  the  role  of 
gravity  in  spin  behavior  and  proposes  scaling  the  models  so  that  the  Froude 
number  (a  dimensionless  quantity  relating  the  relative  effects  of  aero¬ 
dynamic  and  gravitational  forces)  remains  unchanged.  Using  unpowered 
models,  Bihrle  (1976)  also  shows  that  the  type  of  subsequent  spin  motion 
is  highly  sensitive  to  the  aircraft's  initial  condition  (attitude,  control 
setting,  attitude  rates,  etc.);  and  that  spin  motions  which  ensue 
are  highly  sensitive  to  the  sequencing  and  timing  of  the  pro-spin  control 
actions.  He  also  found  that  changes  in  inertias,  side  force  coefficients, 
and  initial  roll  angle  do  not  significantly  influence  developed  spin; 
but  that  the  pitch  damping  coefficient  and  center  of  mass  location  is 
important. 

As  both  Bihrle  (1976)  and  Anglin  (1977)  mention,  it  is  necessary  to 
combine  the  different  types  of  aero  data  in  order  to  have  a  reasonable 
model.  In  most  instances,  rotary  balance  datar  has  limitations  because 
it  is  evaluated  at  relatively  few  control  settings,  and  is  very  un¬ 
reliable  for  angle-of-attack  less  than  about  50°.  In  other  analytic 
results,  Anglin  and  Scher  (1964)  not  only  use  extensively  both  the  con¬ 
ventional  and  rotary  balance  data  to  study  fairly  steady  developed  spins 
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and  recoveries,  but  also  define  and  use  a  non-dimensional  spin-energy 
factor  (Eq.  3.2.48  in  Chapter  III  of  this  report)  to  indicate  the  relative 
difficulty  of  spin  recovery.  This  factor  is  seen  to  be  related  to  the 
antispin  yawing  moment  coefficient.  Further,  they  found  that  the 
antispin  rolling-moment  coefficient  depends  both  on  this  energy  factor 
and  upon  the  moment  of  inertia  about  the  longitudinal  axis. 

Prior  to,  and  concurrently  with,  the  analytic  efforts  briefly  ref¬ 
erenced  above  has  been  an  extensive  program  of  experimental  flight  tests 
and  evaluations  conducted  by  several  government  agencies  and  the  military 
branches.  These  results  are  generally  restricted  to  the  particular  air¬ 
craft  being  studied,  and  are  typically  aimed  at  developing  recovery  tech¬ 
niques  and  avoiding  spin  entry.  Rutan  et  al .  (1970),  Sallada  et  al.  (1967), 
Savidge  (1970),  Glenzer  (1970),  Carlson  (1970),  Krings  and  Weber  (1970) 
and  Shaw  and  Shields  (1970)  all  report  on  the  results  of  spin-oriented 
flight  test  experience. 

Anglin  (1977)  reports  that  much  remains  to  be  done  in  terms  of 
providing  an  aero  data  base  for  spin  regimes  which  will  be  sufficiently 
accurate  to  enable  adequate  simulator  prediction  of  actual  aircraft 
response.  As  we  have  also  found  (Chapter  III),  he  describes  a  large 
region  in  the  yaw  rate  -  angle-of-attack  plane,  located  between  the 
low  angle-of-attack  and  developed  flat  spin  regions,  in  which  neither 
the  conventional  aerodynamics  nor  the  more  recent  rotary-balance  aero¬ 
dynamics  alone  is  adequate  to  describe  the  post-stall  gyrations,  spin 
entry  and  oscillatory  spin  motions  which  characterize  this  region.  Further 
work  is  needed,  he  concludes,  on  understanding  the  behavior  in  this  region; 
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and  this  requires  the  development  of  an  aerodynamic  model  which  incor¬ 
porates  features  of  both  the  conventional  and  rotary  aerodynamics. 

Other  recent  research  worthy  of  note  includes  work  by  Young  (1974), 
who  used  a  steepest  descent  optimization  technique  to  develop  control 
histories  for  spin  recovery.  Adams  (1972)  showed  that  several  spin  modes 
and  types  are  possible  using  representative  aircraft  models,  and  Moore 
et  al .  (1971)  show  that  use  of  a  fixed-base  simulator  can  give  results 
sufficiently  realistic  for  studying  stall/spin  characteristics  of  air¬ 
craft. 

The  previous  work  described  above  has  supplied  us  with  much  of  the 
insight  and  direction  needed  to  adapt  BACTM  to  spin  analysis  problems, 
and  to  clearly  outline  areas  in  which  BACTM  may  be  utilized  to  investigate 
these  problems.  In  the  following  chapters,  particularly  Chapter  III, 
it  will  be  seen  that  our  results  are  in  general  agreement  with  the  above; 
and,  further,  that  BACTM  has  added  fresh  insight  into  many  of  these  problems 
and  possesses  the  capability  to  enhance  even  more  our  understanding  of 
spin  phenomena. 


1.3  Summary  of  Significant  Results 

The  significant  milestones  achieved  on  this  project  are  given  below: 

•  Application  of  BACTM  to  a  six  DOF  aircraft  model,  with  comprehensive 
aero  data  (aircraft  F). 

•  Development  of  simulation  model  for  analyzing  spin  behavior. 

•  Study  of  the  characteristics  of  spin  motion,  and  of  currently  used 
control  procedures  to  simulate  spin  entry  and  effect  spin  recovery. 

•  Representation  of  tabular  aero  data  in  analytical  form  by  means  of 
cubic  and  bi-cubic  spline  functions. 

•  Expansion  and  generalization  of  methodology  for  generating  the 
equilibrium  and  bifurcation  surfaces  of  BACTM.  Reliance  on  para¬ 
metric  continuation  methods  derived  from  work  of  Davidenko  (1953). 

•  Development  of  continuation  methodology,  based  on  work  of  Rhein- 
boldt  (1977)  and  Keller  (1977),  which  can  compute  all  branches 
passing  through  a  "bifurcation  point." 

•  Generation  of  equilibrium  curves  for  aircraft  F  in  flat  spin,  in¬ 
termediate  spin  (angle-of-attack  about  75°)  and  stall  departure 
flight  regimes. 

•  Time  history  simulation  runs  of  the  aircraft  F  model  to 

verify  some  of  the  equilibrium  results;  and  to  begin  an  analysis  of 
the  dynamics  of  spin  entry  and  recovery  from  developed,  flat  spin. 

•  Development  of  an  accurate  and  efficient  means  of  computing  numerical 


derivatives,  using  splines. 
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•  Demonstration  that  flat  developed  spin  motions,  in  which  angle-of- 
attack  is  high  and  pitch  angle  is  low,  may  be  approximated  on  a  first- 
cut  basis  by  assuming  negligible  gravity. 

•  Observation  that  changing  the  value  of  V  by  about  30c;  has  a  more 
significant  impact  on  the  shape  of  the  equilibrium  curves  than  does 
neglecting  gravity. 

•  Observation  that  a  developed  flat  spin  for  aircraft  F  is  featured 
by  an  extremely  tight  spiral,  whose  diameter  decreases  as  the  (pro¬ 
spin)  rudder  setting  gets  more  extreme,  and  which  drifts  slightly, 
due  to  nonsymmetric  lateral  aerodynamics  at  high  a.  The  spin 
velocity,  then,  is  almost  entirely  vertical,  and  the  spin  rotation 
produces  about  0.5  g's  acceleration. 

•  Demonstration  that  a  high-a  limit  cycle  condition  (steep,  oscillatory 
spin)  is  reached  both  by  a  stal 1 -departure  maneuver  starting  from 
trim  conditions,  and  by  a  spin  recovery  maneuver  starting  from  the 
flat  spin  equilibrium  conditions.  The  limit  cycle  family  provides 

the  link  between  the  high-a  equilibrium  states  and  the  trim  equilibrium 
states.  The  recovery  from  this  high  a  limit  cycle  condition  by  a 
fixed  change  in  control  settings  is  possible,  since  the  stability  of 
the  limit  cycle  varies  greatly  over  the  physical  range  of  the  aero- 
surface  control  movements. 
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1.4  Organization  of  the  Report 

The  report  is  organized  into  two  main  chapters,  II  and  III,  and 
a  supporting  chapter,  IV.  Chapter  II  contains  a  discussion  of  the  techniques 
by  which  BACTM  was  modified  to  enable  the  study  of  the  spin  behavior  of 
aircraft;  also  included  in  this  chapter  is  material  on  chaotic  motion 
and  strange  attractors.  Chapter  III  describes  the  use  of  BACTM  to 
analyze  a  particular  aircraft  model.  Aircraft  F,  in  all  high-a  regimes, 
including  spin.  There  is  extensive  discussion  of  Aircraft  F's  behavior 
in  terms  of  equilibrium  surfaces  and  time  history  simulations.  Chapter  IV 
briefly  discusses  other  topics  of  interest,  including  spectral  analysis  of 
chaotic  motions  and  a  preliminary  look  at  using  BACTM  to  synthesize  a 
simple  command  augmentation  system.  Conclusions  anu  recommendations  are 
stated  in  Chapter  V,  and  a  list  of  symbols  and  nomenclature  is  included 
in  Appendix  A. 


CHAPTER  II 


Further  Developments  of  Bifurcation  Analysis  and 
Catastrophe  Thoery  Methodology  (BACTM) 

This  chapter  describes  in  detail  the  numerical  algorithms  used  for 
computing  equilibrium  and  bifurcation  surfaces  for  a  general  nonlinear 
dynamic  model  of  an  aircraft  under  stall  and  spin  conditions.  Notation 
for  the  symbols  presented  in  this  chapter  is  given  in  Appendix  A. 

2 . 1  A  General  Computational  Procedure  Based  on  "Continuation"  Methods 
Continuation  methods  refer  to  those  numerical  techniques  which  "con¬ 
tinue"  a  solution  line,  or  locus,  from  some  point  in  the  state-parameter 
space  where  the  solution  is  known.  That  is,  suppose  the  solution  to  the 
nonlinear  algebraic  system 

f(x,6)  =  0  (2.1.1) 

is  desired.  In  this  equation,  f  and  x  are  each  vectors  of  dimension  n 
(a  more  concise,  mathematical  way  of  saying  this  is  f,xeRn;  which  means 
that  f  and  x  are  elements  of  the  Euclidean  n-space  of  real  numbers,  that 
is,  they  are  n-dimensional ).  Also,  <5  in  this  equation  is  a  scalar,  and 
has  a  special  role  as  the  continuation  parameter.  The  idea  behind  the 
continuation  methods  is  that  we  somehow  know  all  solutions  x  satisfying 
(2.1.1),  for  a  given  <5  =  6q.  These  methods  then  supply  a  means  of  explicitly 
varying  5  from  5q  to  some  desired  value,  6j,  where  analytic  or  numerical 
solutions  are  difficult  to  obtain.  As  an  example,  suppose  we  wish  to  know 
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3  2 

all  values  of  x  satisfying  f(x, 6)  =  x  +2ax  +  bx  +  6  =  0  at  any  given 

value  for  6.  It  is  clear  that  at  6=0,  we  have  three  solutions, 

x=(0,  -a ±  /a2  -  b) .  This,  then,  is  a  natural  starting  point  for  continuing 

the  solution  to  some  nonzero  value  of  6.  Note  in  this  example  that  a  and 

b  are  also  parameters,  but  that  they  remain  fixed  in  value  for  the  entire 

process. 

Continuation  methods  have  been  applied  to  several  varieties  of  prob¬ 
lems  which  are  typically  multi-dimensional,  nonlinear  and  possessing  no 
"analytic"  solution,  i.e.,  a  solution  which  can  be  explicitly  derived. 
Problems  which  have  been  solved  using  continuation  methods  include  certain 
kinds  of  two-point  boundary  value  problems  and  boundary  layer  problems 
(including  singular  perturbation  problems),  in  addition  to  the  algebraic 
problem  defined  by  (2.1.1). 

The  type  of  problem  of  interest  here  is  that  of  solving  a  system  of 
nonlinear  algebraic  equations  of  the  form  (2.1.1).  In  our  applications, 

6  e  {6a,6e,6r} ,  (2.1.2) 

the  set  of  aerosurface  control s--aileron,  elevator,  and  rudder,  respectively. 

Given  that  some  solution  point  x°(6g)  has  been  found  (such  a  solution 
by  definition  satisfies  (2.1.1)),  the  point  x®  is  called  C^-regular  if  the 
Jacobian 


af^x.a) 


(2.1.3) 


an  nxn  matrix,  is  nonsingular  (invertible).  Otherwise,  it  is  a  singular 
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point.  The  notation  C*  refers  to  tie  continuity  of  f(x)  and  all  of  its 
first  derivatives. 

The  solution  to  (2.1.1)  at  x^  is  "continued"  through  other  solution 
points  by  varying  6  over  some  range  of  values.  For  any  regular  point 
(x^,6q)  e  Rn x  R1 ,  the  implicit  function  theorem  ensures  the  existence  of  a 
unique  regular  solution  to  (2.1.1)  through  this  point.  The  notation 
Rn x  R1  means  that  the  (n+l)-space  to  be  considered  is  comprised  of  two 
particular  subspaces:  Rn  for  the  state  variables  x,  and  R1  for  the  scalar 
continuation  parameter.  Continuation  solution  algorithms  qenerally  fall 
into  one  of  two  different  conceptual  classes.  The  first  class  was  initially 
investigated  and  developed  by  Lahaye  (1934,  1948),  and  the  second  approach 
is  usually  attributed  to  Davidenko  (Rail,  1968).  Davidenko's  approach  is 
often  called  continuation-by-differentiation,  and  that  of  Lahaye  belongs 
to  the  class  of  iterative  continuation  techniques.  The  Davidenko  approach 
consists  in  the  application  of  some  suitable  discrete-variable  method  to 
solve 


df  /dx\  9f  n 

ch5  =  FW/  96"  =  5eD,  x(6Q)  =  x  (2.1.4) 

where  D  is  the  set  of  admissible  parameters,  e.g.,  if  6=<5e,  then  D=  [-25,10] 
for  aircraft  F.  Eq.  (2.1.4)  says  that  5  is  to  be  varied  in  a  way  that 
ensures  (2.1.1)  being  always  true.  A  problem  with  this  approach,  es¬ 
pecially  where  n  is  large,  is  the  necessity  for  evaluating  at  each  point 
the  matrix  F  and  the  (nx  l)-vector 


(2.1.5) 
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and  this  is  usually  very  costly.  If  we  try  to  reduce  the  order  of  the 
system  model  to  offset  this  difficulty,  there  may  arise  fatal  accuracy 
problems.  On  the  other  hand,  Lahaye's  iterative  continuation  approach 
uses  a  locally  convergent  iterative  method  (of  a  Newton-Raphson  nature) 
to  solve  (2.1.1),  the  original  equation,  at  an  increasing  (or  decreasing) 
sequence  {6^}  of  parameter  values  in  D.  In  its  basic  form  the  method 
starts  at  the  known  solution  x^1  and  selects  steps  (<$k+1  -  6^)  such 
that  the  last  iterate  at  6k  is  an  acceptable  starting  point  for  the  itera¬ 
tion  at  <5k+1-  At  each  point,  (2.1.1)  is  satisfied  to  within  some  ek  >  0, 
where  e  k  =  ||  f  || .  (See  Appendix  A  for  notation.) 

The  recent  trend  has  been  to  combine  the  two  ideas  by  using  a  feasible, 
low-order  method  of  solving  (2.1.4)  as  a  predictor  and  then  following  it 
with  a  locally  convergent  iterative  process  for  (2.1.1)  as  a  corrector. 

In  particular,  Klopfenstein  (1961),  Keller  (1977),  Rheinboldt  (1977)  and 
Kubicek  (1976)  have  developed  versions  which  seem  particularly  suited  to 
the  task  at  hand:  the  computation  of  equilibrium  and  bifurcation  surfaces 
for  high-performance  aircraft  operating  in  high-a  (nonlinear)  flight  re¬ 
gimes.  One  of  the  major  results  described  in  this  report  is  the  modifica¬ 
tion  and  adaptation  of  the  relevant  techniques  presented  in  the  above 
references  to  two  principal  aspects  of  BACTM  analysis,  the  computation 
of  equilibrium  surfaces  and  bifurcation  surfaces.  The  utilization  of  these 
algorithms,  and  the  refinements  needed  to  handle  certain  situations,  is 
discussed  below. 

More  recently,  we  have  found  that  other  methods  may  be  better  suited  to 

the  particular  application  of  computing  bifurcation  surfaces  for  aircraft 
in  the  spin  regime.  This  is  because  of  both  the  dimensionality  of  this 
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system  (increased  due  to  gravity  coupling)  and  the  highly  nonlinear  charac¬ 
teristics  of  the  motions  in  developed  spin.  Briefly,  these  methods  in¬ 
volve  techniques  which  avoid  the  computation  of  the  Jacobian  (Ralston,  1975, 
Ralston  and  Jennrich,  1978).  A  derivative-free  algorithm  as  well  as  one¬ 
dimensional  search  algorithms  in  the  corrector  phase  of  continuation,  will 
also  be  discussed  below. 

Consistent  with  the  notation  employed  elsewhere  in  this  report,  the 
equilibrium  system  of  equations  has  the  same  form  as  ( 2 . 1 . 1 ) »  viz.: 

f(x,6)  =  0  (2.1.6) 

where  6  e  (6a,<5e,<5r)  and  the  dimension  and  form  of  the  state  x  and  the  function 
(mapping)  f  depend,  in  general,  on  whether  one  is  in  a  spin  or  a  non-spin 
flight  regime.  See  Section  3.2  and  Mehra,  Kessel,  Carroll  (1977),  respec¬ 
tively,  for  the  distinctions.  Eq.  (2.1.6)  is  derived  from  the  aircraft 
dynamic  equations,  which  are  concisely  expressed  as 

x  =  f(x,6)  (2.1.7) 

where 

6  -  (Sa,<5e,6r)  (2.1.8) 

(Hence,  the  mathematical  definition  of  dynamic  equilibrium  is  x  =  0.  For 

dnx 

a  stable  equilibrium,  _ -  =  0,  for  all  n > 0) 

dtn 

The  bifurcation  system  of  equations  has  the  same  form,  but  an  addi¬ 
tional  equation  is  added,  to  specify  the  requirement  that  the  Jacobian 
matrix  F  be  singular  at  a  bifurcation  point.  This  is  in  addition  to  the 
equilibrium  requirement  (2.1.6).  We  denote  the  resulting  set  of  equations  as 
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g(y,6.)  =  o  (2> 

Like  (2.1.6),  (2.1.9)  has  the  form  of  Eq.  (2.1.1),  and  can  thus  be  solved 
by  the  continuation  methods  presented  here.  Furthermore,  (2.1.9)  is  re¬ 
lated  to  (2.1.6)  as  follows: 


where 


and 


where 


5 j » i  ®  (6a,6e,6r) ,  i  f  j; 


S(tf»«f) 


A  £  det(F) 


(2.1.10a) 


(2.1.10b) 


(2.1.1la) 


(2.1.1lb) 


The  solution  of  (2.1.9)  yields  a  curve  6 j ( 6 ^ ) ,  6k  fixed,  in  the  control 
space  called  a  bifurcation  surface.  We  are  at  liberty  to  choose  any  two 
6^  from  the  control  set  (5a,6e,6r),  but  the  third  one,  6^,  remains  fixed 
in  value.  Also,  while  the  bifurcation  surface  is  the  particular  curve 
6^  vs.  Sp  solution  of  (2.1.9)  clearly  supplies  values  for  x  as  well,  at 
each  point  (6,.,Sj)  on  the  bifurcation  surface.  The  system  (2.1.9)  is 
essentially  the  equilibrium  system  (2.1.6)  augmented  by  the  constraint 
A  =  0.  (This  is  discussed  in  more  detail  below.)  The  equilibrium  system 
dimensionality  is  consequently  increased  by  one  in  the  bifurcation  system, 
so  that  an  extra  dependent  variable,  6j e  (6a,6e,6r),  if  i,  may  be  added. 
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The  bifurcation  surfaces  are  consequently  more  difficult  to  generate  for 
two  reasons: 

i)  the  size  of  the  system  is  (n+1),  where  n  is  the  size  of  the 
equilibrium  system.  For  the  non-spin  system  considered  here, 
n=5,  and  for  the  spin  system,  n  =  8.* 
ii)  a  much  worse  problem  than  (i)  is  posed  by  the  presence  of  the 
constraint  on  A,  the  determinant  of  the  Jacobian  of  f.  Even 
in  the  n  =  5  case,  it  is  wholly  impractical  to  expand  A  analyti¬ 
cally,  so  that  evaluating  the  Jacobian  of  g,  (2.1.11a),  of 
necessity  requires  using  a  numerical  differentiation  algorithm 
on  at  least  the  (n+l)th  row  of 


G 


(2.1.12) 


Except  for  extra  core  and  extra  care,  problem  (i)  above  is  relatively 
straightforward  to  surmount.  The  second  problem,  on  the  other  hand,  re¬ 
quires  extreme  caution  and  precision,  in  addition  to  more  complicated 
algorithms.  Consequently,  when  using  continuation-based  methods  re¬ 
quiring  the  first  derivative,  the  core  and  run-time  costs  are  high. 

We  shall  now  discuss  in  more  detail  a  particular  predictor-corrector 
continuation  algorithm  developed  by  Kubicek  (1976). 


2.1.1  Continuation  Method  of  Kubicek  (1976) 

Kubicek' s  method  employs  the  basic  method  of  Davidenko  (Rail,  1968, 


♦There  is  a  dimensionality-flexibility  tradeoff  for  the  spin  system,  which 
centers  around  the  velocity  variable,  V.  V  can  be  solved  for  explicitly 
in  terms  of  the  remaining  7  variables  and  the  controls,  but  at  some  cost 
in  flexibility  and  adaptation  to  several  aircraft  models.  We  have  opted 
more  for  flexibility  in  this  regard,  and  so  n  =  8  for  spin,  not  7. 


22 


Ficken,  1951,  Davidenko,  1953)  to  solve  (2.1.4),  in  combination  with  the 
Newton  method  and  Adams  integration  formulas.  This  particular  method 
has  become  the  basis  for  our  continuation  algorithms,  since  it 
was  found  to  be  capable  of  solving  accurately  and  efficiently  a  wide 
variety  of  nonlinear  algebraic  equations  required  by  the  BACTM  approach 
to  high-performance  aircraft  analysis. 

Kubicek  has  introduced  certain  modifications  to  the  basic  continuation 
methods  of  Davidenko  which  make  the  solution  of  Eq.  (2.1.1)  more  feasible 
on  digital  computers.  In  essence,  this  approach  represents  a  subset  of 
the  methodology  assembled  by  Keller  (1977),  Rheinboldt  (1977),  et  al . ; 
however,  certain  features  of  the  Kubicek  algorithm  are  worth  detailing. 

Basically,  an  arclength  parameter  is  introduced  to  evaluate  the  de¬ 
pendence  x(6),  which  is  assumed  to  be  continuously  smooth  in  the  (n+1)- 
dimensional  space  (x,<5).  (This  assumption  is  not  necessary  in  the  method 
of  Keller,  which  can  handle  the  singularities,  or  bifurcations  points.) 

Quite  often,  x(<5)  is  smooth  in  the  augmented  space  R  1  and  singular  in 
Rn.  In  such  an  instance,  the  system  (2.1.4)  cannot  be  solved,  because  the 
(nxn)  Jacobian  F,  Eq.  (2.1.3),  is  non-invertible,  i.e.,  singular.  This 
happens  at  limit  points,  Fig.  2.1.  At  such  points,  6  is  no  longer  mono- 
tonically  increasing  or  decreasing,  and  F  is  singular.  However,  a  properly- 
selected  arclength  parameter  will  remain  monotonic  at  limit  points;  this 
enables  smooth  continuation  around  limit  points,  as  shown  below. 

By  introducing  the  arclength  parameter  s,  so  that  (x,6)  become 
[x(s),5(s)],  the  system  (2.1.4)  "inflates"  to 

3f 

F5  +  35^  0 


(2.1.13) 


Figure  2. 1 

Bifurcation  and  Limit  Points 
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and  we  choose  initially  a  standard  Euclidean  arclength  relationship 


+  in2  +  62  =  1  (2.1.14) 


where 


%  A  d5 

6  ‘  a? 


(2.1.15) 


Eq.  (2.1.14)  is  comparable  to  the  more  complicated  pseudoarclength  normal i 
zations  introduced  by  Keller  (described  in  the  next  section),  which  are 
useful  in  the  algorithms  which  solve  for  the  branches  at  bifurcation 
points*,  shown  in  Fig.  2.1.  Kubicek  generalizes  the  solution  procedure 
by  generating  a  nonsingular  (nxn)  matrix  of  the  form 


9fl 

3fl 

!!i 

,xk-l'  9xk+l'  '  ' 

3xn+l 

rk  = 


Hjl 

ax 


1 


3Xj  *  * 


af  af 

_ n  _ n 

•  3)w  avr 


af 
_ r 

ax 


n+1 


(2.1.16) 


♦Singular  points  are  points  (x,5)  where  the  Jacobian  matrix  F  is  non-invertible, 
or  singular.  Both  limit  and'bifurcation  points  are  singular  points. 
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In  this  definition,  we  have  set  xn+j=6  for  consistency.  Note  that 

is  nxn  square  because  one  column,  3f/3xk>  is  eliminated.  Since  ke(l,n+l), 

there  are  n+1  possible  Tk  to  analyze. 

We  shall  not  go  into  the  full  detail  here  (Keller,  1977,  has  such 
detail),  but  will  mention  that  at  least  one  nonsingular  rk  does  exist  at 
a  limit  point,  thereby  allowing  continuation  through  that  point.  This 
is  a  consequence  of  the  fact  that,  while  F  is  singular,  it  has  rank  n-1 
at  a  limit  point.  Thus,  the  corank,  equal  to  (n-rank),  is  1.  At  a  bi¬ 
furcation  point,  the  corank  exceeds  1,  and  there  is  no  invertible  rk-  It 
is  possible,  therefore,  to  associate  the  corank  of  F  at  singular  points 
with  the  number  of  equilibrium  branches  intersecting  at  that  point.  At 
a  simple  bifurcation  point,  for  example,  F  has  corank  2  and  two  branches 
intersect  (Fig.  2.1).  At  regular  points,  F  is  nonsingular;  thus  corank 
of  F  is  zero,  but  a  (smooth)  branch  also  passes  through.  If  an  xk, 
l<k<n+l,  can  be  found  so  that  its  corresponding  rk  is  regular,  then 
the  system  (2.1.13)  can  be  recast  as  follows: 


l<j^n+l 


In  this  equation,  the  vector  (dx./ds)  is  n-dimensional ,  as  the  k^1 

J 

element  is  not  used,  and  xn+j-<5.  Eq.  (2.1.17)  can  then  be  solved  for  n 
of  the  n+1  parametric  derivatives: 


l<j$n+l 


(2.1.18) 
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(2.1.20) 


In  Eq.  (2.1.20),  the  sign  ambiguity  is  resolved  by  the  orientation  of  the 
arclength  parameter  s  along  the  curve.  This  may  be  done  somewhat  arbi¬ 
trarily  at  the  solution  starting  point,  (x0>60),  which  must  be  supplied  (or 
obtained  from  a  Newton-type  iteration).  This  solution  branch  will  then 
emanate  in  one  direction  from  the  starting  point,  and  can  be  made  to  emanate 

in  the  other  direction  by  selecting  the  other  sign. 

Thus,  the  method  of  Kubicek  is  more  "robust1’  than  that  of  Klopfen- 

stein  (1961)  in  that  the  latter  retains  the  special  role  of  6  as  the 
continuation  parameter--Klopfenstein  inverts  rn+^,  which  is  actually  F, 
at  every  point  rather  than  the  more  general  T^.  Numerical  difficulties 
are  more  likely  to  be  avoided  by  Kubicek' s  algorithm,  which  at  each  step 
utilizes  the  "best"  continuation  parameter  xk-  The  value  of  the 
subscript  k  is  determined  by  means  of  Gaussian  elimination  using  con¬ 
trolled  pivoting.  That  is,  at  any  point  in  the  reduction  process  for  in¬ 
verting  rk,  the  current  pivot  element  chosen  from  rk,  y.j,  has  the  largest 
magnitude  of  all  candidate  elements.  Once  is  chosen,  all  remaining 
i-h  row  and  j-t-  column  elements  of  Tk  are  eliminated  as  candidate  elements 
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for  selection  of  subsequent  pivot  elements.  Then,  y .  ■  is  used  to  reduce 
all  other  elements  in  column  j  to  zero  by  means  of  the  so-called  "elemen¬ 
tary  matrix  operation" 

{Row  £  of  rk>  =  p£*{Row  i  of  rk)  +  {Row  i  of  rk>  (2.1.21) 

where  the  scalar  is  adjusted  so  that  (2.1.21)  produces  a  zero  for 

element  y^,  if  i.  That  is,  to  zero  all  elements  of  column  j  except  the 

t*  h 

i  — ,  the  right  side  of  (2.1.21)  becomes 
ct'yij  +  vij  ’  °- 

(2.1.22) 

The  process  continues  in  this  manner,  one  column  at  a  time.  Eq.  (2.1.22) 
indicates  the  role  of  the  pivot  elements  y..  •  in  the  matrix  inversion 

process.  If  at  any  step  in  the  process  no  nonzero  y.  .  can  be  found,  the 

1  J 

matrix  is  singular. 

The  process  just  described  is  a  Gaussian  elemination  procedure.  The 
pivoting  is  controlled  in  the  Kubicek  algorithm  by  allowing  each  column 
in  rk  to  be  selectively  scaled.  This  allows  the  user  to  "bias"  the 
selection  of  a  particular  variable  from  (x,6)  as  the  continuation  parameter. 
Typically,  of  course,  the  desired  choice  is  6  =  x  +j.  The  scale  parameter 
for  the  column  associated  with  xn+^  is  then  some  value  less  than  1.,  say 
0.001,  while  those  for  the  other  columns  may  be  kept  at  1.  This  approach 
has  been  developed  by  Deist  and  Sefor  (1967). 

The  complete  procedure  involves  performing  the  elimination  process 
(2.1.21)  n  times  on  the  n-by-(n+l)  array 
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(2.1.23) 


In  so  doing,  n  columns  are  selected  for  r^,  under  the  selection  criterion 
that  y. .  have  the  largest  magnitude  of  the  candidate  elements.  This 
insures  that  the  "most  nonsingular"  is  selected  of  the  (n+1)  possi¬ 
bilities.  The  one  unselected  column  in  this  process  becomes  column  k. 
When  the  resulting  reduced  r  matrix  is  deprived  of  this  column,  re¬ 
sults.  rk  is  of  (2.1.16),  but  operated  on  several  times  by  the  ele¬ 
mentary  matrix  operations  (2.1.21).  ?k  has  the  property  that  there  is 
exactly  one  nonzero  element  (indicated  by  x)  in  each  of  its  n  rows  and 
columns,  e. g. , 


(2.1.24) 


Eq.  (2.1.17)  thus  becomes 


rk*  +  \  xk  =  0 


(2.1.25) 


where  f  is  the  k^1  column  of  the  reduced  r  matrix. 
xk 

Upon  expansion,  (2.1.25)  has  the  form  of  (2.1.19),  so  that  the  Bj 
from  the  latter  expression  are  readily  extracted,  and  the  n  x^*dx./ds 
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are  found  from  (2.1.19),  i  f  k.  Also,  (2.1.20)  produces  x^.  Given  the 
n+1  x  elements,  the  Newton  step  size  Ax  is  obtained  for  the  full  (n+1)- 
dimension  x.  Note  that  Ax^ =  0  in  this  procedure.  What  the  above  process 
has  effectively  done  is  compute  a  Newton  iteration  step,  based  on  the 
standard  multi -dimensional  Newton-Raphson  formula 

A*(P)  =  -r-1f(x^),  (2.1.26) 

where  p  is  the  iteration  counter  (x^  is  given),  and  x  includes  the 
original  x,  plus  5. 

To  sumnarize,  we  have  solved  (2.1.17)  for  the  x.-dx./ds,  j/k ,  l<j<n+l; 
i.e.,  we  have  found  the  Bj  in  (2.1.19).  Using  Eq.  (2.1.20),  which  represents 
the  standard  Euclidean  arclength  relationship  utilized  by  Klopfenstein  and 
Kubicek,  all  (n+1)  quantities  (x,6)  are  determined.  These  derivatives  are 
then  used  by  the  Kubicek  algorithm  to  predict  the  next  point  on  the  curve. 
This  is  done  by  using  Adams-Bashforth  integration  formulae.  The  logic  of 
Kubicek  and  Deist  and  Sefor  which  regards  the  parameter  6  as  interchangeable 
with  any  state  element  x.  at  each  point  makes  this  algorithm  very  robust 
in  terms  of  singularities  and  numerical  roundoff  difficulties.  This  is 
because  the  value  k  can  change  from  point  to  point  along  the  continuation 
solution.  After  the  Adams  integration  step  (predictor  step),  the  new 
point  (x,6),  where  all  n+1  quantities  have  changed,  is  the  next  starting 
point  for  the  Newton  iteration  (corrector  step)  to  the  point  again  satis¬ 
fying  (2.1.6). 

The  Kubicek  algorithm  has  been  discussed  in  considerable  detail  be¬ 
cause  it  plays  a  major  role  in  the  adaptation  of  BACTM  to  handle  large 


order,  complex  systems.  However,  we  have  implemented  other  methods  which, 
although  not  used  as  much  at  this  point,  do  exhibit  promise  in  certain 
aspects.  Ultimately,  it  is  hoped  to  develop  a  unified,  comprehensive  and 
flexible  package  for  BACTM  which  utilizes  the  most  appropriate  algorithm 
for  the  task  at  hand.  Some  of  the  other  algorithms  which  are  currently 
being  developed  will  now  be  discussed. 


2.1.2  Continuation  Method  of  Davidenko  (1953)  and  Numerical  Differentiation 


If  Eq.  (2.1.1)  can  be  put  in  the  form: 


dx 

^  =  M(x,6)  (2.1.27) 


any  number  of  numeric  integration  methods  will  be  able  to  solve  for  a 
complete  branch  of  solutions  once  an  initial  solution  is  found.  This 
transformation  is  accomplished  by  differentiating  Eq.  (2.1.1)  to  obtain 
(2.1.4)  so  that 


When  this  equation  can  be  solved,  either  an  interpolator  integration 
method  (such  as  Runge-Kutta)  or  a  predictor-corrector  method  (such  as 
Euler-Newton  or  Adams-Bashforth)  can  be  used  to  continue  along  a  branch. 
This  is  the  essence  of  Davidenko' s  method,  which  has  subsequently  been 
refined  considerably. 

The  first  task,  then,  is  to  solve  for  the  partial  derivative  of  f 
with  respect  to  any  x^  or  6,  since  these  comprise  elements  of  F^.  In 
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the  past  we  have  solved  for  these  derivatives  analytically  while  setting 
up  the  problem.  However,  this  severely  limits  the  range  of  problems  which 
may  be  solved  and,  moreover,  makes  the  process  of  applying  the  technique 
to  a  new  problem  extremely  complicated. 

As  a  result,  we  have  implemented  a  numeric  differentiation  routine. 

This  routine  is  based  on  a  cubic  spline  fit*  to  the  function  f(x,6), 
evaluated  at  selected  points  y  =  (x,6),  centered  on  yQ,  the  point  where  the 
actual  derivative,  3f/3y  =  3f/3(x,6)  -  T  is  desired.  To  understand  the 
process  more  readily,  consider  the  scalar  case:  x,  f,  6eR'.  The  goal, 
then,  is  to  numerically  evaluate  3f/3x  and  3f/36  at  yQ=  (xg,dg).  For 
3f/3x,  evaluate  f  at  the  five  points 

f((xQ  +  pe),60),  p  =  -2, -1,0, 1,2  (2.1.29) 

Note  that  we  have  f(xQ,6Q)  at  p =  0.  The  increment  size  e  is  chosen  so  that 
llf(tf0>-f(tf)ll  *  10"4  (2-1-30) 

This  choice  of  e  allows  sufficient  accuracy  of  the  fit  without  introducing 
serious  numerical  difficulty,  e  must  provide  a  large  enough  spread  in  Ay 
so  that  the  slope  obtained  is  representative,  yet  not  be  so  large  as  to 
deteriorate  precision.  We  then  use  each  of  the  five  y  as  knots  for  the 
spline  fit.  The  polynomial  approximation  resulting  from  this  is  analytic 
at  yg,  by  definition,  so  that  evaluation  of  3f/3x  at  yg  is  straightforward. 

Similarly,  for  3f/35  at  yg,  fix  x  at  Xg,  evaluate  f  at  the  five  points 

f(x0,(6g+pe))  (2.1.31) 


♦Curve-fitting  using  splines  is  discussed  in  detail  in  Section  3.2.1. 
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and  make  a  spline  fit  and  evaluation  as  before. 

The  extension  to  x,f e  Rn  follows  immediately.  For  each  9f/3xi , 

make  five  evaluations  of  the  vector  f  by  varying  x.  only,  keeping  x-  =  x.  , 

1  J  J0 

j^i,  and  6  =  6g.  Each  of  the  fm  are  then  spl  ine-  fitted,  and  this  pro¬ 
duces  the  set  of  vectors  (3f  /3x,),  for  each  x.,  and  6. 

mi  i 

2.1.3  Keller-Klopfenstein  Continuation  (Keller,  1977) 

Limit  points.  Using  the  Davidenko  method,  all  points  along  a  solution 
branch  can  be  computed  as  long  as  the  nxn  matrix  F  can  be  inverted.  If 
F  is  singular  with: 

dim  N(F)  >  1  (2.1.32) 

where  N(*)  denotes  null  space,  special  procedures  must  be  used  at  such  a 
point.  The  null  space  dimension  here  is  equivalent  to  the  corank  of  F. 
Adding  a  "normalization"  equation  and  a  new  parameter  to  the  system  can 
avoid  this  problem  in  some  cases.  Keller  (1977)  uses  the  normalization: 

N3  $  9xJ(x-x0)  +  (1-0)SO(«-«O)  -  (s-sQ)  -  0  (2.1.33) 

where  s  is  the  arclength  parameter,  0  <  0  <  1 ,  and  Xq,  <$q,  and  Sp  are  the 
values  at  the  initial  point  for  the  branch  in  question.  The  system  now 
becomes: 


(2.1.34) 


where : 
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Y  = 


(2.1.35) 


This  system  can  be  dealt  with  in  the  same  manner  as  Eq.  (2.1.1). 

The  normalization  is  said  to  have  better  numerical  properties  near 
bifurcation  points.  Keller  proves  that,  using  (2.1.33)  as  the  normali¬ 
zation  equation,  3g/3y  is  nonsingular  if  and  only  if: 

F  is  nonsingular  (2.1.36) 


or: 


3f/36  t  R(F),  where  R(*)  denotes  range  space.  (2.1.37) 

Case  (2.1.37)  corresponds  to  a  "limit  point."  At  such  a  point  (see  Fig¬ 
ure  2.1)  two  branches  do  not  intersect,  but  dx/d6^«>.  In  this  case,  solu¬ 
tion  of  the  augmented  system  (2.1.34)  continues  normally.  The  proof  can 
be  developed  by  utilizing  Gaussian  elimination  techniques.  We  shall  dis¬ 
cuss  this  and  other  aspects  of  Keller's  method  in  more  detail  in  later 
reports. 

For  dealing  with  limit  points  alone,  Keller's  method  is  somewhat 
cumbersome,  particularly  in  the  choice  of  for  an  arclength  relationship. 
The  method  of  Kubicek,  discussed  in  Section  2.1.1,  is  quite  adequate  at 
limit  points,  using  the  simpler  normalization  (2.1.14).  However,  the 
Keller  approach  is  more  comprehensive,  and  can  handle  the  computation  of 
equilibrium  branches  at  bifurcation  points.  We  shall  now  outline  how  this 
is  done,  saving  some  of  the  detail  for  later  reports. 

Bifurcation  points.  If,  at  some  point  on  a  branch,  neither  condition 


_*****■>■> 
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(2.1.36)  or  (2.1.37)  is  satisfied,  3g/3y  will  be  singular  and  continuation 
of  the  augmented  system  will  also  fail.  This  is  a  "bifurcation"  point, 
where  at  least  two  branches  cross.  At  such  a  point  two  steps  must  be 
accomplished.  First,  the  bifurcation  must  be  skipped  over  so  that  con¬ 
tinuation  may  proceed  along  the  same  branch.  Second,  a  point  on  the  other 
branch  must  be  found  as  an  initial  point  for  continuation  along  it. 

Any  predictor-corrector  method  can  be  used  to  skip  over  a  bifurcation 
point.  This  simply  involves  choosing  an  initial  solution  and  step  size 
for  a  single  prediction  step  for  which  the  correction  step  will  converge 
to  a  solution  past  the  bifurcation  point. 

Finding  the  second  branch  at  the  bifurcation  point  is  more  complicated. 
To  begin,  consider  a  simple  bifurcation  point,  at  which  two  branches  inter¬ 
sect.  Here,  the  rank  of  T  is  n-1.  r  is  an  n-by-(n+l)  matrix,  recall. 
Multiple  bifurcation  points  have  more  than  two  branches  intersecting  at  y*, 
and  the  rank  of  r  is  less  than  n-1.  By  examining  the  nature  of  the 
various  terms  in  the  power  series  expansion  of  f(x,<5)  near  such  a  bifur¬ 
cation  point,  y*,  one  realizes  that  the  two  branches  emanating  from  y*  lie 
tangent  at  y*  to  a  plane  defined  by  the  two  eigenvectors  associated  with 
the  two  zero  eigenvalues  of  rTr,  which  is  a  square  (n+l)-by-(n+l)  matrix 
of  rank  n-1.  This  matrix,  then,  has  a  corank  of  (n+l)-(n-l)  =  2;  hence, 
two  zero  eigenvalues.  Once  r  is  found,  therefore,  the  plane  is  readily 
determined,  and  all  branches  can  be  located  by  a  search  in  this  plane 
for  f=0  points.  This  search  is  done  along  an  arc  of  fixed  radius  from 
y*,  and  sufficiently  close  to  y*. 

To  summarize,  the  algorithm  of  Keller  differs  from  that  of  Kubicek 


in  these  significant  areas: 
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(i)  The  pseudo -arc length  normalization  of  Keller  is  more  compre¬ 
hensive  than  the  purely  geometric  relationship  (2.1.14)  employed  by 
Kubicek,  and  is  capable  of  dealing  with  bifurcation  points; 

(ii)  The  Gaussian  elimination  procedure  for  selecting  as  con¬ 
tinuation  parameter  at  each  point  seems  to  be  capable  of  effecting  numerical 
solutions  more  efficiently.  Incorporating  it  into  Keller's  algorithm 
could  well  enhance  its  usefulness. 

(iii)  A  final  distinction  which  has  practical,  if  not  theoretical, 
significance,  is  the  means  by  which  one  proceeds  along  the  solution 
branches.  Both  schemes  utilize  predictor-corrector  algorithms,  with  some 
form  of  Newton  method  as  a  corrector: 


-new  -old 


(2.1.38) 


However,  Kubicek  uses  the  Adams-Bashforth  explicit  multistep  method,  with 
variable  order,  to  integrate  (2.1.19)  and  (2.1.20),  whereas  Keller  suggests 
a  modified  Euler-Newton  scheme,  which  is  really  a  first-order  Adams 
algorithm.  There  is  more  flexibility  in  the  more  complete  Adams  method. 

2.1.4  Applications  of  the  Kubicek  Continuation  Algorithm 

The  algorithm  works  exceptionally  well  as  coded  for  solving  aircraft  F 
equilibrium  surfaces,  (2.1.6),  both  in  spin  and  non-spin  regimes.  There 
is  consequently  no  need  to  discuss  the  algorithm  per  se  with  regard  to 
these  surfaces.  However,  the  computation  of  bifurcation  surfaces  of  air¬ 
craft  F,  defined  by  (2.1.9),  did  necessitate  more  care  in  setting  up  the 
problem  and,  in  the  case  of  spin  bifurcation  surfaces,  required  in  addition 
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some  modification  to  the  algorithm  itself.  These  points  will  now  be 
discussed  in  more  detail. 

2. 1.4.1  Bifurcation  Surfaces 

Let  x(6)  represent  an  equilibrium  solution  to  (2.1.6).  It  is  obvious 
that  the  variation  of  one  element  of  6  (6=  (6a,6e,6r))  will  generate  an 
equilibrium  surface  in  the  state-control  space.  A  bifurcation  point  on 
the  equilibrium  surface  implies  a  change  in  structural  stability  for 
control  variations  in  the  neighborhood  of  the  bifurcation  point.  Points  A 
and  B  are  bifurcation  points  in  Figure  2.2,  and  the  loci  of  their  pro¬ 
jection  onto  the  control  subspace  is  what  we  call  a  bifurcation  locus. 


Figure  2. 2 
Equilibrium  Surface 


The  set  of  bifurcation  loci  on  a  particular  control  subspace  is  called  a 
bifurcation  surface.  The  locus  is  generated  by  varying  any  two  of  the 
three  controls,  holding  the  third  one  fixed.  It  is  obvious  that  the  bi¬ 
furcation  points  are  a  subset  of  the  equilibrium  points.  The  criterion 
for  an  elementary  bifurcation  point*  is 

*An  elementary  bifurcation  point  has  a  zero  eigenvalue.  More  general  bi¬ 
furcation  points  such  as  Hopf  Bifurcation  points  have  purely  imaginary 
eigenvalues.  (See  Mehra  et  al_.  (1977)  for  details.) 
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A  ^  det(  F)  =  det 


0 


(2.1.39) 


(i.e.,  the  Jacobian  is  singular).  The  difficulty  at  bifurcation  points 
arises  from  the  fact  that,  if 


f  ( x ,  6 )  =  0 


(2.1.6) 


is  true,  then 

dx  3f 

Fd?  0  <2-' 

when  6  is  one  of  the  controls  selected  as  a  parameter.  (See  Section  2.1.) 
It  is  seen  from  Fig.  2.2  that  dx/d6  is  the  slope  of  the  f=0  locus  for 
values  of  x  and  6  on  that  locus.  Also,  points  A  and  B  are  characterized 
by  the  fact  that  dx/d<5,  the  slope,  is  infinite.  Hence,  the  continuation 
solution 


dx 

d6 


(2.1.40) 


breaks  down.  This  is  equivalent  to  saying  that  the  inverse  of  the 
Jacobian  F  does  not  exist,  i.e.,  Eq.  (2.1.39)  holds. 

Thus,  unadulterated  continuation  methods  break  down  at  such  points, 
as  these  methods  solve  for  x(<5)  by  integration  of  (2.1. 40).  Kubicek's 
algorithm  avoids  this  by  introducing  an  arclength  parameter  and  by  aug¬ 
menting  the  Jacobian  with  an  extra  column  representing  the  parameter,  and 
eliminating  (via  Gaussian  reduction)  that  column  which  leaves  the  most  nonsin¬ 
gular  square  matrix  (this  amounts  to  interchanging  the  parameter  6.  with 
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an  x^  as  necessary;  xk  then  becomes  the  parameter  and  continuation  via  an 
equation  like  (2.1.40)  or  (2.1.18)  remains  valid). 

To  summarize,  Eqs.  (2.1.6),  evaluated  for  some  starting  value  x^,  and 
(2.1.40),  which  continues  the  solution  from  Xq,  generate  equilibrium 
surfaces.  Bifurcation  surfaces--represented  by  points  a  and  b  in  the  con¬ 
trol-space  in  Figure  2.2  --are  generated  in  a  similar  way,  with  the  basic 
equilibrium  system  being  enhanced  by  one  dimension  (representing  the  con¬ 
straint  (2.1.39)).  Thus,  the  bifurcation  system  is  indeed  given  by  (2.1.9). 

For  equilibrium  surfaces,  one  of  the  controls  is  selected  as  a  parameter, 
leaving  the  other  two  fixed.  For  bifurcation  surfaces,  one  of  the  controls  is 
still  a  parameter,  but  one  of  the  remaining  two  controls  becomes  a  state  variabl 
because  of  the  introduction  of  a  in  g  (see  (2.1.11a)).  The  bifurcation 
surfaces  then  are  generated  by  a  system  similar  to  (2.1.13),  viz.: 

.  d9  ngnty  /ag  \ds 

9  "  d?  '[lyTJds  +  ^as'.Jds-  '  0 

or,  in  compact  form, 

9  “  +  96  *  =  0 

As  outlined  in  Sec.  2.1.1,  Eq.  (2.1.42)  is  solved  for  the  (n+1)  derivatives* 

£  (the  (n+1)^1  element  of  ^  is  5.,  as  defined  in  (2.1.10a)),  as  functions 
of  the  scalar  6.  This  latter  derivative  is  then  determined  from  the  arc- 
length  normalization  relationship 

*In  the  spin  flight  regime,  n  =  8.  Thus,  the  equilibrium  system  (2.1.6)  is 
8th  order,  the  bifurcation  system  (2.1.9)  is  9th  order,  and  the  matrix  r 
of  Sec.  2.1.2  is  a  (9x10)  array. 


(2.1.41) 


(2.1.42) 
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ft  +  C&)2  =  1  (2.1.43) 

See  Sec.  2.1.1  for  details. 

2. 1.4. 2  Non-spin  bifurcation  considerations 

In  the  non-spin  case,  n  =  5  (gravity  effects  are  neglected  and  V  is 
assumed  constant).  However,  even  though  this  system  is  considerably  smaller 
than  the  spin  system  (for  which  n  =  8),  the  r  array  is  of  size  6x7.  Fur¬ 
thermore,  because 

A  =  det(F)  (2.1.11b) 


represents  the  (n+1)  —  element  of  g  (hence  row  (n+1)  of 


G  = 


3yj 


), 


(2.1.12) 


there  are  serious  computational  problems  to  consider. 

These  problems  center  on  the  computation  of  G  and  the  (n+l)-by-(n+2)  array 


(2.1.44) 


which  is  often  more  difficult  than  inverting  the  submatrix  of  r.r^.  For 
example,  the  bottom  row  of  r,  (2.1.44),  is  given  by 


(2.1.45) 


where  y  is  defined  by  (2.1.10a).  It  is  of  practical  necessity  to  com¬ 
pute  at  least  this  row  using  a  numerical  differentiation  algorithm. 
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Even  for  n =  5,  deriving  the  analytic  expression  for  a  is  excessively  labor¬ 
ious,  not  to  mention  all  of  its  partials;  needless  to  say,  there  is  no 
point  in  speculating  about  the  n  =  8  (spin)  case. 

Thus,  the  last  row  of  r  in  either  the  spin  or  the  non-spin  case  is 
determined  by  numerical  differentiation.  The  first  column  element  in 
this  row,  Y(n+i)  i  is  given  by  3A/3p.  As  described  in  Sec.  2.1.2,  the 
numerical  differentiation  of  this  term  involves  five  evaluations  of  A(p), 
with  all  other  variables  fixed.  For  each  state,  and  the  two  selected 
controls,  this  must  be  done;  altogether  seven  times  for  each  of  the  last 
row  elements  of  T,  in  the  non-spin  case.  Thus,  35  evaluations  of  A  are 
needed  each  time  an  evaluation  of  r  is  made.  There  can  be  several 
evaluations  of  r  made  for  each  point  on  the  continuation  solution,  due 
to  the  iterative  nature  of  the  Newton-corrector  steps.  Every  evaluation 
of  A  requires  full  evaluation  of  the  matrix  F.  It  is  possible  to  do  this 
using  numerical  differentiation,  but  there  is  obviously  a  tremendous 
saving  in  time  to  be  had  if  the  terms  in  F  can  be  analytically  derived, 
as  well  as  all  other  elements  of  r  (F  is  a  submatrix  of  T)  for  which 
this  is  feasible. 

Thus,  the  above  strategy  was  adopted  for  computing  r,  both  in  the  spin 
and  non-spin  cases;  i.e.,  use  analytic  expressions  for  the  elements 
as~  far  as  possible,  using  numerical  differentiation  only  for  the  last 
row  of  r,  (2.1.45).  This  modification  provided  the  opportunity  to  compare 
the  precision  of  the  numerical  differentiation  results  with  the  "exact" 
expressions,  and  the  numerical  adequacy  of  the  former  was  verified. 

Another  modification  made  to  run  bifurcation  solutions  was  to  evaluate 


the  aerodynamic  coefficients  only  once  for  every  point  actually  accepted 
as  a  solution  point,  and  not  every  time  F,  or  even  f,  is  evaluated.  This 
results  in  considerable  time  savings,  as  the  aero  data  for  aircraft  F 
requires  interpolation  routines. 

The  modifications  discussed  above,  when  applied  to  the  non-spin  bi¬ 
furcation  system,  can  generate  almost  a  solution  point  per  CPU  second  on 
the  CDC  6400,  an  improvement  by  about  a  factor  of  50  on  the  unmodified 
system. 

2. 1.4. 3  Spin  Bifurcation  Surfaces 

The  spin  bifurcation  system  worsens  the  "curse  of  dimensionality." 

In  this  regard,  incorporating  an  algebraic  soluti  r  for  V,  which  reduces 
the  basic  system  dimension  from  8  to  7,  does  not  help  very  much  computa¬ 
tionally.  This  is  because  all  of  the  partials  involving  V  in  F  would 
have  to  be  carried  along  by  the  chain  differer  iation  rule.  V  could  be 
specified  to  be  constant,  but  it  is  as  yet  unclear  whether  this  action 
would  obscure  transition  dynamics  in  the  pos'  -stall  departure  and  spin 
entry  flight  regimes. 

With  this  background,  the  n  =  8  spin  system  was  incorporated  into  the 
BACTM  bifurcation  package,  modified  as  described  above  for  the  smaller, 
non-spin  system.  Numerical  difficulties  were  encountered  immediately, 
which  centered  around  the  Newton-Raphson  iteration  algorithm.*  More 
specifically,  the  problem  lies  with  the  Newton  corrector  step,  computed 
by  the  relation 


A*(k)  -  (*(k+1)  -*(k))  - 


(2.1.46) 


♦These  difficulties  were  present  as  well  when  a  scaled  velocity,  n-  v/vscaie 


was  introduced.  V 


scale 


was  given  a  value  which  put  n  in  the  range 


of  values  of  the  other  variables,  all  of  which  are  in  radian  units. 
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where 


l  =  (p.q,r,a,B,V,6,((),6j)T  =  (x,6^)T 


S  -  9(^’6i) 


(T’) 


(2.1.10a) 

(2.1.11a) 


and 


A  -  det 


=  det  (F) 


(2.1.11b) 


(In  the  continuation  process  employed  by  BACTM,  the  Newton-Raphson  al¬ 
gorithm  itself  is  considered  a  Corrector,  as  it  refines  the  approximation 
g-0  iteratively.  The  Prediction  step  is  made  by  the  Adams-Bashforth 
multi  step  method.)  The  algorithm  above  works  quite  well  for  all  of  the 
BACTM  systems  except  the  spin  bifurcation  system.  In  this  case,  if 
is  defined  as 

£(k)  -  lla(Y(k),<5i)||  (2.1.47) 

where  y^  is  the  k^1  Newton  iterate  at  a  (predicted)  point  designated  by 
( 6^ , 6 j ) ;  then,  it  typically  happens  that 

£(k+l)  >  £(k)  (2.1.48) 

for  the  spin  system.  Worse,  this  happens  over  a  range  of  k  values,  so 

(kl 

that  the  sequence  {Ay'  '}  becomes  "unstable"  with  respect  to  the  itera¬ 
tion.  If  a  solution  g  =  0  is  found  under  these  conditions.,  it  invariably 


is  at  a  point  far  removed  from  the  regime  of  interest. 

By  the  nature  of  the  Newton  algorithm,  assuming  ge  C1,  the  class  of 

functions  continuous  in  all  first  derivatives,  there  exists  a  region  R 
m 

centered  on  y'  '  such  that 

£*  <  (2.1.49) 

where  the  stepsize  producing  e*  relates  to  Ay^  by 

Ay*  =  p(k^Ay(k)  (2.1.50) 

O 

(  k) 

where  p'  ;e(0,l]  is  a  Euclidean  metric  between  two  points  in  R.  Basically, 

this  says  that  if  one  cuts  back  sufficiently  on  the  step  size  as  gi''en  by 

Eq.  (2.1.46),  then  one  will  eventually  find  a  step  size  (2.1.50)  for  which 

(2.1.49)  holds,  if  g  is  "sufficiently  smooth."* 

The  above  represents  an  aspect  of  one  general  approach  to  the  basic 

f  kl 

goal  of  ensuring  a  decreasing  sequence  {e'  that  is,  the  one-dimensional 
search  algorithm.  These  methods  accept  the  minimizing  direction  as  com¬ 
puted  by  (2.1.46),  even  though  it  uses  only  first  order  information. 

Then,  a  scalar  one-dimensional  search  is  made  along  this  direction  for  a 

point,  y*,  which  is  "minimizing."  Several  algorithms  exist,  and  two 

f  kl 

have  been  implemented.  The  first  of  these  merely  halves  Ay'  ' ,  evaluates 

(kl 

g  at  that  point,  compares  its  norm  to  e'  and  halves  the  stepsize 

a^ain,  continuing  until  either  e*<e^  or  1 1 Ay j |  <  10 The  latter  con¬ 
fix 

dition  indicates  either  that  g  is  not  continuous  at  y'  ' ,  or  that  the 

first  order  information  at  that  point  is  inadequate.  This  algorithm 

was  successful  in  some  regions,  but  less  so  in  others,  even  with 

*This  is  easy  to  see  geometrically  if  one  constructs  a  scalar  system  from 
Eq.  (2.1.46). 
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modifications  such  as  permitting  (2.1.48)  up  to  three  consecutive  times. 

Another  one-dimensional  algorithm  which  was  implemented  involves 

f  kl 

considering  only  the  segment  along  the  direction  Ayv  ’  between  the  points 
(k)  (k)  (k) 

'  and  '+Ayv  If  one  defines  a  scalar  arclength  parameter  a 
(kl  fk+ll 

such  that  y(a =  0)  =  y'  '  andy(a?l)=yv  ,  then  it  is  possible  to  con¬ 
struct  a  cubic  polynomial  in  a  using  only  previously  computed  quantities: 
g(y^)),  G(y^)  and  g(y^+^).  The  cubic  is  then  assumed  to  approximate 
the  function  e  in  this  interval,  and  its  minimum  value  should  be  close 
to  the  function  minimum.  This  method  has  fared  little  better  than  the 
first  one  mentioned,  but  is  has  not  been  tested  completely.  A  combination 
of  the  above  two  algorithms  is  also  being  considered. 

Other  algorithms  under  consideration,  but  as  yet  not  implemented,  in¬ 
clude  one-dimensional  search  and  a  least  squares  method  based  on  step¬ 
wise  regression.  The  latter  is  useful  especially  if  G  is  nearly  singular. 
If  singularity  problems  arise,  the  method  of  stepwise  regression  may  be 
employed,  as  it  eliminates  variables  so  that  g  is  better  parameterized. 

The  remaining  variables  are  those  which  "do  the  most  good"  in  reducing 
2  2 

e,  or  e  =  ||g||  .  Stepwise  regression  can  deal  with  nonlinear  formula¬ 
tions  and  bounded  spaces. 

Another  general  approach  which  shows  great  promise,  but  is  yet  to 
be  tested,  involves  derivative-free  algorithms.  In  these,  the  Newton 
approach  is  forsaken  entirely  in  favor  of  methods  which  rely  only  on  the 
evaluation  of  g{y).  By  the  proper  sequence  of  (n+1)  such  evaluations, 
a  secant  plane  in  n-space  may  be  constructed.  At  each  iteration,  the  use 
of  this  plane  along  with  the  "optimality  criterion"--!' .e. ,  minimizing 
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the  distance  between  a  point  in  this  plane  and  the  desired  point--provides 
information  for  replacing  one  of  the  (n+1)  y-points  with  a  new  point. 

Then,  the  new  set  of  (n+1)  points  defines  a  new  secant  plane  which  im¬ 
proves  the  optimality  criterion.  In  a  convergent  sequence,  the  secant 
plane  iterates  in  the  limit  to  the  tangent  plant  which  contacts  the  function 
g  at  y*,  such  that  g(y*)  =  0,  the  desired  solution.  Ralston  (1975)  and 
Ralston  and  Jennrich  (1978)  discuss  a  specific  derivative-free  algorithm 
(D.U.D.--doesn' t  use  derivatives)  which  holds  some  promise  for  our 
applications. 
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CHAPTER  III 

BACTM  Applied  to  the  Study  of  Aircraft  F  Stall -Spin  Behavior 

This  chapter  describes  in  detail  the  aircraft  model,  designated 
Aircraft  F,  which  is  used  in  this  report  for  the  study  of  high-a,  post¬ 
stall  and  spin  motions.  Also  discussed  in  detail  are  descriptions  of 
stall  and  spin  behavior;  the  nature  of  spin  and  stall;  the  dynamic  equa¬ 
tions  which  describe  the  above  motions;  the  aerodynamic  data  for  aircraft 
F;  and  BACTM  results  for  aircraft  F  in  spin  entry,  stalling  maneuvers, 
wing  rock,  post-stall  gyrations,  developed  spin  motion,  and  spin  recovery. 
An  explanation  of  the  spin  behavior  of  aircraft  F  completes  the  chapter. 

3 . 1  Nature  of  Aircraft  Stall  and  Spin  Behavior 

In  this  section  a  brief  overview  will  be  presented  to  provide  some 
insight  into  the  physical  phenomena  of  aircraft  stall  and  spin.  Briefly, 
an  aircraft  encounters  a  stall  condition  when  loss  of  lift  occurs  due  to 
excessive  buildup  of  angle-of-attack  (a).  High  performance  aircraft  are 
particularly  susceptible  to  this  phenomenon  if  only  because  their  design 
goal  is  to  operate  near  the  extremes  of  the  flight  envelope  in  accom¬ 
plishing  mission  objectives.  Stall  and  spin  are  related  because  post¬ 
stall  behavior  can  include  departure  into  spin  conditions.  Since,  as 
will  be  shown  in  this  report,  spin  motions  can  be  stable  equilibria,  it 
is  important  to  understand  enough  about  spin  phenomena  to  be  able  to 
effect  recovery  control  sequences  which  transfer  the  aircraft  state  from 
the  domain  of  attraction  of  spin  equilibrium  points  into  the  domain  of 
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attraction  of  nonspin  stable  equilibrium  points. 

3.1.1  Nature  and  Characteristics  of  Spin 

For  most  aircraft,  high  angle-of-attack  departures  from  stall  quite 
often  proceed  to  entry  into  one  of  several  spin  modes.  Whether  a  spin 
condition  is  achieved,  and  if  so,  which  one,  depends  to  a  great  extent 
on  the  particular  aircraft  configuration  and  control  settings,  as  well  as 
the  flight  condition.  In  general,  the  spin  modes  are  characterized  by 
the  incidence  angle  (angle  of  attack),  a.  A  typical  classification  of 
"erect"  spin  modes  has  five  categories  (Rutan,  et  al .  (1970)): 

(1)  Steep  -  Smooth 

(2)  Steep  -  Mild  Oscillation 

(3)  Steep  -  Oscillatory 

(4)  High-a  -  Highly  oscillatory 

(5)  Flat 

For  certain  configurations,  "inverted"  spin  is  possible,  but  this 
particular  phenomenon  has  not  been  investigated  in  this  report.  Other 
investigators,  e.g.,  Adams  (1972)  and  Young  (1974),  classify  erect  spin  as 
steep,  intermediate,  or  flat,  and  oscillatory  or  steady.  At  any  rate, 
the  higher  the  equilibrium  angle  of  attack  is,  the  flatter  the  spin. 
Typical  values  of  a  during  a  fully  developed  spin,  derived  from  a  study 
of  F-4  behavior  (Adams  1972),  are  a  =  80° -85°  for  a  flat  spin,  =72°, 73° 
for  an  intermediate  spin,  and  =45°-  60°  for  a  steep  spin. 

There  are  basic  characteristics  of  spin  motion  common  to  all  of  the 
spin  modes  defined  aoove: 

(1)  The  angle  of  attack  remains  greater  than  the  stall  angle 
of  attack  (aSTALL); 
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(ii)  the  center-of-mass  follows  a  helical  path,  with  the  net 
velocity  being  almost  totally  vertical  and  constant  sink  rate; 

(iii)  the  aircraft  attitude  changes  at  a  steady  rate  (of  magni¬ 
tude  w),  and  the  axis  of  rotation  is  almost  totally  vertical. 

In  the  development  of  the  equations  for  analyzing  developed,  or  equilibrium, 
spin,  the  above  characteristics  will  be  incorporated  as  dynamic  constraints. 
Other  assumptions,  such  as  constant  speed  (V),  will  be  discussed  later. 

The  visual  cue  for  a  spin  is  excessive  and  continuous  yaw  rate. 

The  fully-developed,  or  equilibrium,  spin  is  often  achieved  after 
only  a  few  rotations.  Once  the  spin  is  established,  the  trajectory  is 
essentially  vertical.  In  this  situation,  equilibrium  is  a  result  of  a 
balance  of  the  aerodynamic  (decomposed  into  lift  and  drag),  gravity  and 
the  centrifugal  forces  (the  latter  arising  from  the  helical  motion). 

The  drag  vector  opposes  the  velocity,  and  so  is  largely  vertical;  hence, 
the  lift  vector  tends  to  lie  in  the  horizontal  plane,  and  the  radius  of 
the  helix,  R,  adjusts  until  the  resultant  centrifugal  term,  w2R,  balances 
the  lift  term.  During  an  established  spin,  the  presence  of  non-zero 
sideslip,  6,  is  quite  often  prominent  in  the  generation  of  the  coupling 
moments  which  maintain  the  equilibrium  spin.  These  particular  moments 
are  usually  affected  by  the  rudder  setting. 

When  an  aircraft  enters  a  spin  from  a  basically  straight  and  level 
flight  condition,  the  center  of  gravity  follows  a  path  that  initially 
was  horizontal,  but  changes  to  a  vertical  spiral.  Such  a  significant 
change  in  flight  condition,  caused  basically  by  entry  into  a  stall  region, 
is  bound  to  produce  at  least  transient  oscillatory  behavior  in  p,q,r,  the 
rotation  rates  about  the  vehicle  roll,  pitch  and  yaw  axes,  respectively. 
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Generally,  such  oscillations  dissipate  before  one  or  two  spin  revolutions 
have  occurred;  however,  in  some  instances,  these  oscillations  may  grow  in 
amplitude.  Characteristic  of  such  oscillatory  spins  is  the  large  change 
of  aircraft  attitude  with  respect  to  the  horizon  as  roll  rate  p  oscillates 
(Kerr (1956)).  The  oscillations  may  cause  the  motion  to  change  from  a  spin 
to  a  post-stafl  gyration*,*  in  which  rolling  motion  is  prominent. 

We  have  devoted  the  bulk  of  our  spin  analysis  efforts  in  this  report 
to  the  study  of  spin  behavior  when  the  aircraft  is  in  flat  spin  situa¬ 
tions.  This  is  primarily  because  flat  spins  have  tended  to  be  the  most 
troublesome  ones,  in  that  recovery  from  them  is  usually  very  difficult 
to  achieve.  This  is  because  much  of  the  spin  equilibrium  regime  is 
stable,  which  requires  active  control  for  recovery.  Aircraft  susceptible 
to  post-stall  entry  into  flat  spins  tend  to  have  such  characteristics  as 
lengthened  fuselage  forebodies,  increased  relative  distribution  of  mass 
in  the  fuselage,  and  wings  with  short  span.  The  design  ideal  is  that  the 
aircraft  equilibrium  be  unstable  in  flat  spin  over  as  wide  a  control  regime 
as  possible,  so  that  recovery  from  flat  spin  is  easy. 

3.1.2  Nature  and  Characteristics  of  Stall 

Stall  is  the  condition  of  dramatic  loss  of  lift  due  to  a  change  in 
the  operating  state  of  the  aircraft.  Stall  entry  is  typically  a  longi¬ 
tudinal  phenomenon,  in  that  application  of  elevator  will  cause  angle- 
of -attack  to  grow  excessively.  Lift  is  proportional  to  angle-of-attack 


until  a  value  for  a  is  reached. 


“STALL’  at  which  ^ow  seParati°n  abound  the 


wings  occurs.  The  wake  becomes  turbulent,  lifting  capacity  is  sharply 


reduced,  and  the  aircraft  is  said  to  have  stalled.  Stall,  then,  is  a 


major  phenomenon  of  high-a  flight.  Often,  pre-stall  buffeting  or  wing  rock 
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warn  the  pilot,  though  not  always  in  time.  Subsequent  to  stall,  depending 
on  the  state  of  the  aircraft,  there  may  be  post-stall  departure  into  spin, 
roll  departure,  or  autorotation;  all  of  these  are  undesirable  motions, 
particularly  spin,  which  can  be  an  equilibrium  state,  structurally  stable 
over  the  available  range  of  control  values. 

The  pilot  generally  has  to  assume  active  control  of  the  vehicle  once 
the  aircraft  has  stalled,  and  there  may  be  extreme  conditions  which  only  an 
autopilot  can  deal  with.  In  any  event,  one  of  the  main  goals  of  aircraft 
and  flight  control  system  design  is  to  extend  the  operational  flight  regime 
as  much  as  possible  by  avoiding  stall  situations,  and  to  provide  the 
pilot  with  adequate  warning  when  he  is  about  to  encounter  stall. 
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3 . 2  Spin  Equations  of  Motion,  Assumptions,  Constraints 

The  equations  of  motion  which  adequately  model  spin  behavior  require 
both  the  full  set  of  coupled  nonlinear  terms  and  a  data  base  for  the 
aerodynamic  coefficients  which  is  comprehensive  enough  to  include  regimes 
in  which  spins  occur.  In  the  equations  that  follow,  engine  thrust  and 
gyroscopic  terms  are  not  included.  The  role  of  engine  thrust  in  spin 
entry  and  recovery  dynamics  will  be  studied  in  subsequent  projects.  Pre¬ 
vious  work  (Grafton  1966,  Lusby  1961)  has  shown  these  effects  to  be  small;* 
furthermore,  the  greatly  reduced  airflow  along  the  longitudinal  axis  may 
cause  serious  engine  damage,  erratic  thrust  behavior,  and  flameouts.  Finally, 
the  variation  of  atmospheric  density  p  with  altitude  h  is  neglected,  so 
that  dynamic  pressure  q  is  a  function  only  of  airspeed  V.  The  basic  equa¬ 
tions,  then,  are: 


+  \-(&c  -a 


;ine)s 


a  =  9  +  I  - 1  ^  cx  -  y-  sine  +  r  sinejsina 
+  ("mV’  cz  +  V  cose  C0Si^  "  P  sine^cosaj 


sees 


cosa 


6  =  -[01  Cx  -  V  sine)s1n6  +  >] 

+  (ml  Cy  +  V  cos0s1n<^)cos^ 

-  j(|y-  Cz  +  ^  cos8cos<t>)sinB  -  pjsi 


na 


(3.2.1) 


(3.2.2) 


(|)  =  (If  Cx  '  V  sin9)cosacos0  +  (jjjf  cy  +  v  cosesin<t»)sinB 
+  (m?  cz  +  f  cosecos^sinacose 


(3.2.3) 


♦However,  thrust  seems  to  play  a  prominent  role  in  spin  entry  dynamics. 
See  Sec.  3.3. 
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(3.2.4) 


(3.2.5) 


(3.2.6) 


These  equations  are  derived  principally  from  Adams  (1972),  but  some 
of  the  terms  in  the  expansion  of  the  force-moment  coefficients  were  added 
from  other  sources  (Moore,  Angl in  (L97l)  and  Brady(l969)).  Notation  is  pre¬ 
sented  in  Appendix  A. 

The  following  kinematical  relationships  are  needed  to  fully  describe 
the  motion: 


e  =  q  cost|>  -  r  sin<f> 

(3.2.7) 

4>  =  p  +  q  tanesin<}>  +  r  tanecos<t> 

(3.2.8) 

^  =  a  sin<j)sece  +  r  cos^sece 

(3.2.9) 

R  =  V^cosn  +  V^sinn 
a  y 

(3.2.10) 

Rn  =  -vjsinn  +  V^cosrt 
x  y 

(3.2.11) 

(3.2.12) 
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where 

V*  =  V[cosa  cose  cos<|*  cose  +  sina(cos4i  sine  sin^  -  sin^  cos4> ) 

+  sina  cos6(cosii<  sine  cos$  +  si n<^  sin<j>)]  (3.2.13a) 

Vy  =  V[cosa  cose  sim|»  cose  +  sina(simp  sine  simji  +  cos^  cos <j>) 

+  sina  cose(sin^  sine  cos$  -  cos^  s i n<f> ) ]  (3.2.13b) 

-•  •  »  •  .  •  •'  - 

V*  =  V[-cosa  cose  sine  +  sina  cose  sin$  +  sina  cose  cose  cos<f>]  (3.2.13c) 


In  Eqs.  (3.2.1)  to  (3.2.6)  the  aerodynamic  coefficients  are  expanded 
as  follows: 


Cx  "  Cx(a,B,6  =  0)  +  Cx  (a,B)6e 

6e 


(3.2.14a) 


Cy  =  C^(a,e,5  =  0)  +  cu  (a,e)6e  +  Cu  (a,e)6a  +  Cv/  (a,e)6r 


r6e 


6a 


6r 


(zv)[Cy  (“ 


)p  +  Cw  (a)r 


(3.2.14b) 


Cz  =  Cz(a,e,4=0)  +  Cz  (a,6)6e 

6e 


(3.2.14c) 


-  (^(a, 6,6  =  0)  +  (a,e)<$e  +  CA  (a,e)6a  +  (a,e)6r 


6e 


6a 


6r 


h(w)  k  (a)p +  s (a)r 


(3.2. 14d) 


Cm  =  Cm(a‘6,6  =  0)  +  Cm^(a,8)6e  +^wjC^(a)q 


(3.2. 14e) 
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C 


n 


C  (a,6,6  =  0)  +  c  (a,6)6e  +  C  (a,B)<5a  +  C  (a,6)6r 
<5e  n6a  n6r 


Cn  (a)P  +  C  (a)r 
P  r 


(3.2. 14f) 


The  above  system,  Eqs.  (3.2.1)  -  (3.2.12),  allows  a  complete  time  his¬ 
tory  for  general  motion  of  the  aircraft.  Note  that  Eqs.  (3.2.1)  -  (3.2.8) 
are  a  self-contained  sub-system,  not  requiring  information  from  Eqs.  (3.2.9)  - 
(3.2.12);  the  converse  is  not  true. 

We  shall  now  digress  a  bit  to  discuss  how  the  equilibrium  system  used 
by  BACTM  is  developed  from  the  dynamic  system  of  equations.  By  "equi¬ 
librium"  is  meant  dynamic  equilibrium.  That  is,  we  are  seeking  all  sets 
of  points  (x,6)  for  which  f(x,6)  =  0.  Here,  the  elements  of  the  vector  f 
consist  of  the  right-hand  terms  in  Eqs.  (3.2.1)-  (3.2.8),  for  one  system; 
the  elements  of  the  vector  x  are  (p,q,r,a,B,V,0  ,<j>) ,  for  that  system;  and 
the  elements  of  <5  are  (<5a,<5e,<5r) .  At  f=0,  of  course,  all  time  derivatives 
of  x  are  zero,  so  that  p,q,r,  etc.,  have  fixed  values.  This  is  dynamic 
equilibrium.  There  remains  the  very  important  issue  of  classifying  the 
equil ibrium  points.  First,  some  definitions:  A  stable  eguiliprium  point 
is  one  in  which  motions  originating  at  some  point  in  the  neighborhood  of 
the  equilibrium  point  (x,6),  ultimately  return  to  y;  an  unstable  equi¬ 
librium  point  is  one  in  which  any  motion  beginning  at  a  non-equil iDrium 
point  near  y  will  ultimately  diverge  from  All  equilibria  for  linear 
systems,  x=Fx+G<5,  are  either  stable  or  unstable.  For  nonlinear  systems, 
such  as  we  are  dealing  with  here,  motions  not  only  may  either  converge  or 
diverge  asymptotically,  but  also  may  develop  into  1 imit  cycles.  This  is 
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a  motion  in  which  the  time  solution  x(6;t),  also  referred  to  as  an  orbit 
or  trajectory,  returns  at  some  time  (t  +  T)  to  its  value  at  t,  where  T  is 
fixed.  Depending  on  its  amplitude,  frequency,  and  other  issues  (some  of 
them  probably  subjective),  a  limit  cycle  may  or  may  not  be  considered  a 
stable  motion.  There  are  other  points  to  emphasize,  however.  First, 
limit  cycles  themselves  are  stable  or  unstable*,  in  the  same  sense  that 
equilibrium  points,  or  curves,  are;  that  is,  motions  starting  near  a 
stable  limit  cycle  ultimately  end  on  it,  and  conversely  for  unstable 
limit  cycles.  Second,  a  point  on  a  stable  limit  cycle  is  not  an  equi¬ 
librium  point,  because  x^O. 

Equilibrium  points  are  classified  typically  by  investigating  the 

first  order  term*  in  the  power  series  (e.g.,  Taylor)  expansion  of  f(x,6), 

evaluated  at  y.  By  the  nature  of  nonlinear  systems,  then,  the  validity 

of  the  stability  classification  is  restricted  to  those  regions  centered 

on  y  where  the  zeroth  and  first  order  terms  of  the  expansion  of  f  are 

almost  equal  to  f  itself.  (This  comment  has  particular  relevance  to  limit 

cycle  motions.)  Specifically,  the  eigenvalues  of  the  Jacobian  matrix,  F  or 

3f(x,6)/9x,  evaluated  at  y=  (x,6),  are  used  to  classify  the  equilibrium 

point.  For  the  eighth-order  dynamic  system  (3.2.1)  -  (3.2.8),  F  is  8x8, 

and  yields  eight  eigenvalues.  If  all  of  these  have  negative  real  parts 

(any  complex  eigenvalue  has  a  conjugate  mate,  as  all  coefficients  in 

the  equation  x =  Fx  are  real)  at  y,  that  equilibrium  point  is  stable;  if 

one  or  more  has  a  positive  real  part,  the  equilibrium  is  unstable.  As 

alluded  to  above,  this  classification,  depending  as  it  does  only  on  linear 

analysis,  cannot  provide  quantitative  information  about  limit  cycles 

*It  can  be  proven  that  first  order  analysis  is  sufficient  for  accurate 
classification  of  "nice"  nonlinear  functions  if  done  in  a  region  close 
to  the  point  about  which  the  expansion  is  done. 
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(we  are  developing  for  BACTM,  however,  a  novel  approach  for  quantitative 
analysis  of  limit  cycles,  using  continuation  methods.  This  will  be  de¬ 
scribed  in  more  detail  later).  In  a  region  where  limit  cycle  behavior 
exists,  the  "governing"  equilibrium  branch  will  be  unstable,  typically 
possessing  one  very  lightly  damped  complex  pair  of  eigenvalues  (i.e., 
real  part  close  to  zero,  but  positive);  the  others  have  negative  real 
parts.  Thus,  when  the  (linear)  eigenanalysi s  shows  a  lightly  damped 
complex  pair,  all  we  can  say  is  that  limit  cycle  motion  is  expected  in 
the  region.  We  cannot  quantify  the  limit  cycle  (e.g.,  amplitude,  period, 
stable  or  unstable,  etc.)  without  further  analysis.  The  equil ibrium 
point  which  indicates  limit  cycle  motion  is  actually  an  unstable  one, 
and  this  is  in  fact  true  in  a  "local"  sense  (that  is,  in  the  region 
about  y  where  the  linear  approximation  to  f  is  valid).  Motions  starting 
near  such  a  point  diverge;  however,  once  their  amplitudes  are  large 
enough  so  that  nonlinear  influences  are  greater,  a  limit  cycle  may 
result.  If  this  happens,  the  global  motion  is  stable,  and  is  the  asymp¬ 
totic  limit  of  motions  emanating  near  the  locally  unstable  equilibrium 
point.  Such  a  limit  cycle  is  a  stable  attractor,  and  we  have  seen  many 
examples  of  these  (see  Mehra  et  al_. ,  1977).  Unstable  limit  cycles  do 
play  an  important  role  as  well,  but  their  quantification  is  considerably 
more  difficult  than  that  of  stable  limit  cycles.  Global  bifurcations 
deal  with  the  annihilation  of  stable  limit  cycles  by  unstable  ones,  as 
system  parameters  (e.g.,  control  settings)  change. 

In  summary,  then,  equilibrium  points  may  be  stable  or  unstable.  For 
a  certain  type  of  unstable  point,  the  one  with  one  unstable  complex  pair 
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of  eigenvalues,  limit  cycle  notion  may  result,  but  not  locally.  Motion 
is  locally  divergent  from  all  unstable  points  and  locally  convergent  to 
all  stable  points. 

We  return  now  to  the  aircraft  spin  equilibrium  system.  In  order 
to  find  the  spin  equilibrium  points,  two  extra  conditions  are  needed: 
i)  an  equilibrium  is  specified  a  priori  by  the  requirement 


P 


q 


r 

a 


=  0 


8 


V 


This  result  is  derived  from  the  equilibrium  requirement  for  constant 
angular  velocity  and  steady  helical  motion. 

ii)  the  spin  characteristics  are  specified  as  dynamical  constraints, 
to  be  incorporated  in  the  system  Eqs.  (3.2. 1)  -  (3.2. 12) .  These  constraints 
are: 


h=h*,  a  specified  constant;  (3.2.15) 

this  decouples  the  h  equation  and  lets  q=q(V)  only. 

a?  =  5/z1  (3.2.16) 

where  <*>  represents  the  total  angular  velocity  of  the  aircraft; 


ip(t*)  =  0 


(3.2.17) 
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where  t*  is  the  time  at  which  the  equilibrium  solution  is  made.  This 
relationship  results  in  no  loss  of  generality,  because  of  the  natural 
decoupling  of  ip  from  the  basic  dynamical  system,  Eqs.  (3.2.1)-  (3.2.8). 

Final ly , 

V  =  R^T  -  hz1  (3.2.18) 

where  T  is  the  unit  vector  in  the  direction  tangent  to  the  trajectory,  and 
z*  is  the  locally  vertical  unit  vector.  This  is  the  constraint  for  helical 
motion.  It  leads  to  the  relationships 

V  •  ft  =  0 

Vn  T  2  T  2 

R  =  -r~  .  V„  =  r  +  r  (3.2.19) 

15-1  H  x  y 

•  • 

n  = 

In  (3.2.19),  is  the  magnitude  of  the  horizontal  component  of  velocity, 
and  ijj  is  the  heading  rate.  When  constraints  (3.2.15)-  (3.2.17)  are  in¬ 
corporated  into  the  full  dynamic  system  of  equations,  the  system  reduces 
to  a  five-dimensional  set  of  nonlinear  algebraic  equations  for  the 
equilibrium  points  in  the  state  space.  Before  presenting  this  equilibrium 
system,  a  few  consequences  of  incorporating  the  constraints  will  be  de¬ 
tailed.* 

The  requirement  =  0  allows  for  the  direct  algebraic  solution  of  V  as 
a  function  of  (a,B,6,p,r,e,<J>) ,  using  Eq.  (3.2.3).  Specifically,  V  is  the 
solution  of  a  quadratic  equation.  Condition  (3.2.16)  results  in  the 
identities 

♦Later,  we  will  show  an  eighth  order  spin  equilibrium  which  was  actually  used 
for  the  numerical  solutions.  This  system  was  used  because  it  is  less  con¬ 
figuration-dependent  and  because  the  continuation  algorithms  developed  for 
BACTM  (Ch.  II)  can  efficiently  handle  large  order  systems. 
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4>  =  0 
0  =  0 
ii  =  w 

(3.2.20) 

p  =  -oi  sine 
q  =  u  cose  sin<j> 
r  =  o)  cose  cos<t> 

Considering  the  last  three  of  these  identities,  it  is  seen  that  a  new 
state  variable,  w,  has  been  defined,  replacing  the  set  (p,q,r),  for  a  net 
reduction  of  two  in  the  order  of  the  system.  Also  recall  that  the  time 
rate  of  change  of  the  Euler  angles  (\p , 0 ,d>)  is  not  an  orthogonal  set. 

This  is  not  true  of  the  roll,  pitch,  and  yaw  rates  (p,q,r),  which  are 
orthogonal  body  axis  components  of  the  total  aircraft  angular  velocity 
vector,  oj.  As  vector  components,  ('l»,e,|)  do  add  up  to  uj  also,  but  only  in 
special  cases  may  be  associated  with  yaw  rate,  0  with  pitch  rate,  and 
so  on.  This  means,  as  (3.2.20)  shows,  that  all  Euler  angle  rates  but  ip 
may  be  zero,  yet  because  of  projections,  all  of  (p,q,r)  are  nonzero. 

The  reader  needing  further  amplification  of  this  aspect  of  aircraft 
kinematics  is  referred  to  an  appropriate  text  covering  kinematics,  e.g. , 
Goldstein  (1950),  and  in  particular,  Etkin  (1972). 

The  equilibrium  system  of  equations  is: 

0  =  a  =  w  s i n<j>cose  +  [-sina(F^Cx  -  ^  sine  +  u  cosecos^sine) 

+  cosa(F^Cz  +  ^  cosecos<()  +  w  sinesine)]sece 


(3.2.21) 
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B  -  -[sina(F1Cx  -  2.  sine)  +  w  cosecos$]cosa 

+  cose^Cy  +  ^  cosesin<j>)  -  [sine(F1Cz  +  |  cosecos$)  +  u>  sine]sina 


=  (cosecosijiu)  = 


"M/(' 


I.  -  I 


g  — *-5 -  sinecosesin<f> 


-  /l  +  -Zj~  ^  ^  Fg  w2cos2esin4>cos4> 


(3.2.22) 


+  b  F4^5C£  +  Cr 


(3.2.23) 


(-sinew)  =  ~~  -oj2  cosesin<j>cosij> 
9  . 


I  -  I 


I  -  I 


+  w2  Fg  h  -  JL. 


+  b  F2<C*  +  W 


F8  +  ^T 


sinecosesin<j> 


(3.2.24) 


0  =  (cosesinifiw)  =  c  Fg  Cm  +  u)2jF^(cos2ecos2(j>  -  sin2e) 


sinecosecos4> 


(3.2.25) 


The  system  of  Eqs.  (3.2.17)  -  (3.2.21)  may  be  compactly  expressed  as 


f(x,6)  =  0 


(3.2.26) 


where 


jj  -  (ct, 8,u),0 ,<t> ) 


(3.2.27) 


In  the  above  system,  the  following  identities  are  introduced 
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F  =  evs. 

rl  2m 


f  - 

h2  “  21. 


F  -  £X£S 
F3  2Iy 

F  =e^l 

h4  21 


F  =  J« 
5  I 


_  xz 


F3  =  F5F6 


F  =  — 
10  2V 

F  = 
rll  2V 


(3.2.28) 


Note  the  presence  of  V  in  the  F. .  This  requires,  then,  the  introduction  of 
the  quadratic  expression  for  V  into  the  system,  which  is  usually  no  prob¬ 
lem  for  solution  procedures  which  are  iterative  in  nature,  given  a  starting 
point.  When  the  system  is  solved  for  the  equilibrium  values  (a,ij,w,0,£) , 
and  V,  then  other  quantities,  such  as  R,  may  be  found: 


R  = 


V 

w 

0) 


[sina(cos$  -  cosesin^)] 


It  is  worthwhile  also  to  monitor,  at  the  equilibrium  solution  points, 
the  following  parameters,  which  have  been  found  by  experimenters  to  be  use¬ 
ful  in  detecting  departure  into  spin: 

Directional  Departure  Parameter 

C  =  C  COSa  -  C  sina  (3.2.29) 

Vn  6  x  6 


Lateral  Control  Directional  Parameter 


LCDP  =  C  -  C. 

n3  *6 


Cn  +  K  Cm 

n6a  6r 


+  K  C„ 


6a 


' fir 


(3.2.30) 


62 


where  K  =  — -  . 

Both  C  and  LCDP  should  be  positive  for  stable  motion.  However, 
^DYN 

they  were  originally  intended  to  compensate  for  the  dynamic,  or  forced 

oscillation,  derivatives.  Therefore,  when  more  complete  aero  data 

bases  are  utilized,  parameters  such  as  C  and  LCDP  may  lose  some 

eDYN 

significance.  They  are,  in  essence,  a  preliminary  design  tool.  Other 

derivatives  such  as  C  and  C  may  be  interesting,  and  it  is  certainly 

n3 

of  interest  to  verify  that 
a  >  aSTALL‘ 

Because  of  the  improved  methodology  developed  during  the  project, 
it  has  become  possible  to  use  a  somewhat  more  general  set  of  dynamic 
equations  for  the  study  of  spin  behavior.  This  set  is  an  eighth  order 
system,  and  it  has  the  advantage  of  flexibility  and  ease  in  terms  of  de¬ 
riving  the  elements  of  the  system  Jacobian  matrix  analytically.  This  system 
consists  of  Eq.  (3.2.1)  -  (3.2.6)  as  well  as  the  kinematic  relations  (3.2.7) 
and  (3.2.8),  and  has  been  coded  directly.  By  reasonable  choice  of  initial 
conditions,  the  solution  of  this  set  either  for  time  history  trajectories 
or  for  equilibrium  surfaces  will  automatically  incorporate  the  constraint 
relationships  (3.2.15)  -  (3.2.17). 

Additionally,  a  second  set  of  velocity  state  variables  is  used  when 
time  history  solutions  are  generated.  This  set  uses  the  body-axis  com¬ 
ponents  of  V,  namely  (u,v,w),  as  state  variables  instead  of  the  set  (V,a,e) 
which  is  defined  by  Eqs.  (3.2.1)  -  (3.2.3).  The  set  (V,a,B)  is  the  wind- 
axis  velocity  state  and  the  set  (u,v,w)  is  the  body-axis  velocity  state. 


63 


whose  dynamical  relations  are  given  by: 

u  =  -g  sine  +  vr  -  wq  +  ^-  Cx  (3.2.31) 

v  =  g  cosesin<(!  +  wp  -  ur  +  ^-  (3.2.32) 

w  =  g  cos9cos<j)  +  uq  -  vp  +  ^-  Cz  (3.2.33) 

where  C  ,  C  ,  C  are  the  total  aero  force  coefficients  along  the  aircraft 
x  y  z 

x- ,  y-  and  z-axes.  An  auxiliary  set  of  relations  enables  computation  of 


variables  of  importance: 

a  =  tan_1(w/u) 

(3.2.34) 

6  =  sin_1(v/V) 

(3.2.35) 

V  =  /u2  +  v2  +  w 2 

(3.2.36) 

The  use  of  the  body-axis  set  is  disadvantageous  in  that  the  aero¬ 
dynamic  coefficients  are  functions  of  a  and  6,  so  that  (3.2.34)  and  (3.2.35) 
must  be  carried  along.  However,  with  this  set,  and  with  the  equation  for 
heading  angle  (also  the  yaw  angle  in  a  yaw-pitch-roll  inertial -to-body 
Euler  transformation  sequence)—  Eq.  (3.2.9)— several  kinematic  and  dynamic 
variables  of  interest  in  the  force-moment  equations  may  be  generated 
rather  easily.  Also,  this  set  is  "cleaner"  in  form  than  the  wind-axis 
set,  as  it  avoids  the  transcendental  terms  in  a  and  6  which  represent  the 
body-to-wind  axis  transformation. 

Some  more  terms  which  are  worth  monitoring  in  the  study  of  stall 


and  spin  behavior  include 
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TURNS  = 

Ztt 


E  =  h  mV2  +  H  to* I*u 


E  =  mVV  +  I  pp  +  I  qq  +  I  rr 
x  y  z 


where 

V  =  (uu  +  vv  +  ww)/V 
if  the  body-axis  velocity  set  is  used 

a)  =  /p2  +  q2  +  r2 

VVERT  =  ~u  sin0  +  ^  sin<Jl  +  w  C°Sc|) )COS0 , 

positive  downward, 

VNORTH  =  cos'i;tu  cose  +  (v  sin4»  +  w  cos<j>)sine] 
-  sini|<(v  cos4>  -  w  sin<(>) 

^EAST  =  cose  +  (v  s'*n<i)  +  w  cos<j>)sine] 

+  cosjp(v  cos<p  -  w  sin<*>) 


Es  =  Js  I/2/(qSb) 


(3.2.37) 

(3.2.38) 

(3.2.39) 

(3.2.40) 

(3.2.41) 

(3.2.42) 

(3.2.43) 

(3.2.44) 

(3.2.45) 

(3.2.46) 

(3.2.47) 

(3.2.48) 
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where 


Iy  =  sin2elx  +  (sin2^  +  cos24>Iz)cos2e  +  Ixzcos4>sin2e  (3.2.49) 


(when  Ixy  =  lyz  - 

“spin 


(3.2.50) 


E  is  a  scalar  quantity  representing  total  vehicle  kinetic  energy 
(assuming  a  purely  rigid  body);  the  first  term  is  due  to  translational 
(center  of  mass)  motion,  and  the  second  represents  the  rotational  contri¬ 
bution.  Here,  I  is  the  moment  of  inertia  tensor,  written  as  a  square 
3x3  matrix;  thus,  pre-  and  post-multiplication  of  I  by  the  angular 
veloc.cy  vector  oj=  (p,q,r)  results  in  the  scalar  quantity  (assuming  that 
the  off-diagonal  terms  of  I,  I  ,  I  ,  and  I  ,  are  zero) 

"  Xa  j  c 

*SU>«I-W  =  hlj>2  +  hlv q2  +  hi  r2  (3.2.51) 

Note  that  the  only  sensible  axis  system  for  coordinatization  of  I  is  the 

Sody-axis  system,  and  thus  w  is  also  in  this  system  for  compatibility. 

E  is  merely  the  time  rate  of  change  of  E,  the  kinetic  energy.  Dspi N 

is  a  Euclidean  metric  which  is  an  indicator  of  how  far  a  given  point 

in  the  state  space,  x,  is  from  a  known  (input)  stable  equilibrium  spin 

location,  In  (3.2.50),  x^  is  a  component  of  the  state  vector  x, 

say  a,  x  is  the  value  of  that  component  at  a  known  spin  equilibrium 
si 

condition,  and  n  is  the  number  of  state  elements  in  x.  DsP IN  *  then*  is 
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intended  to  provide  some  indication  of  how  far  a  particular  point 
x e  Rn  is  from  a  known  spin  equilibrium  point,  x£. 

Equations  (3.2.37)  -  (3.2.50)  represent  relations  for  quantities 
which  are  of  interest  in  investigating  spin  motion,  as  well  as  the  quan¬ 
tities  defined  by  Eq.  (3.2.29),  (3.2.30)  and  C  .  Also,  selected  terms 
from  the  basic  dynamic  set  of  equations  are  usually  of  interest. 

3.2.1  Aircraft  F  Configuration;  Representing  Aero  Data  by  Spline  Functions 

Because  of  the  completeness  and  compactness  of  data  which  is  relevant 
to  the  study  of  spin  conditions,  aircraft  Configuration  B  from  Moore, 

Anglin  (1971)  will  be  used  first  in  this  study.  We  shall  henceforth  call 
this  configuration  aircraft  F.  Additionally,  there  are  several  sources 
of  aerodynamic  data  for  the  F-4  Series  of  aircraft,  namely,  Rutan  (1970), 

Adams  (1972),  Moore  and  Anglin  (1971),  Brady  (1969);  however,  some  extra 
effort  is  required  to  coalesce  and  reduce  these  data  to  the  form  of  air¬ 
craft  F,  and  so  the  F-4  Series  will  be  investigated  later.  Moore,  Anglin 
(1971)  use  the  F-4  data,  but  not  for  the  study  of  equilibrium  spins,  so 
no  high-a  data  is  supplied  there.  Data  for  Aircraft  F  are  presented  in 
Tables  I  and  II. 

Aircraft  F  is  a  variable  sweep  fighter  aircraft  whose  aerodynamic  data 
is  somewhat  equivalent  to  that  of  the  F— 111 ,  although  it  must  be  emphasized 
that  values  were  modified,  especially  the  Cn  derivatives,  to  allow  the  simu¬ 
lation  results  to  readily  produce  "typical"  spin  motions. 

It  has  been  decided  to  model  the  aero  data  for  aircraft  F  by  using 
cubic  and  bi-cubic  splines.  A  cubic  spline  is  a  means  of  representing  data 
points  by  third  order  polynomials.  The  values  of  the  four  coefficients 
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which  specify  the  polynomials  are  determined  by  continuity  and  smoothness 
conditions  at  the  so-called  "knots,"  which  are  typically  but  not  neces¬ 
sarily  the  data  points  themselves.  Thus,  the  cubic  polynomial  between 
any  two  knots  differs  from  that  between  the  neighboring  knots,  but  their 
respective  coefficients  are  selected  so  that,  at  their  common  knot,  the 
value  itself  and  (for  cubic  splines)  the  value  of  the  first  two  deriva¬ 
tives  match.  The  bi-cubic  spline  is  conceptually  similar  to  the  cubic 
spline,  except  that  it  is  a  cubic  polynomial  in  two  variables,  and  not 
one.  See  Eq.  (3.2.52).  We  have  developed  analytic  functions  in  (a, 6) 
for  the  coefficients  by  using  cubic  and  bi-cubic  splines,  for  the  follow¬ 
ing  reasons: 

(i)  splines  assure  smoothness,  especially  at  the  boundaries  (knots); 

(ii)  The  numerical  algorithm  employed  by  BACTM  to  generate  the  equi¬ 
librium  surfaces  requires  the  partial  derivatives  of  the  right 
hand  side  of  Eqs.  (3.2.1)-  (3.2.6)  with  respect  to  the  state 
and  control  variables--i .e. ,  it  is  required  to  generate  the 
matrix* 

*3f  ‘ 

.3*  . 

f"B) 

Since  (m,e)ex,  expressions  such  as  — r— 1 —  will  be  required. 

“  dot 

The  stability  analysis  in  the  neighborhood  of  the  equilibrium 
surfaces  also  requires  this  matrix. 

(iii)  The  data  need  not  be  supplied  over  uniform  increments.  This 
has  the  potential  for  saving  much  core  on  the  computer,  since 
the  smoother  portions  of  the  data  don't  require  as  many  points 
to  adequately  define  the  function. 

*The  coefficients  are  linear  in  the  control  variables  6  so  this  aspect  is 
straightforward,  since  the  partials  are  the  control  stability  derivatives. 
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TABLE  I 

Aircraft  F: 

Weight,  N  (lb) .  222  410  (50  000) 

Wing  Area,  m2  (ft2) .  48.77  (525.0) 

Wing  span,  m  (ft) .  19.20  (63.0) 

Mean  aerodynamic  chord,  m  (ft) .  2.76  (9.04) 

lx,  kg-m2  (slug-ft2) .  67  790  (50  000) 

ly,  kg-m2  (slug-ft2) .  427  348  (315  200) 

lz,  kg-m2  (slug-ft2) .  476  564  (351  500) 

Ixz,  kg-m2  (slug-ft2) .  0  (0) 

Maximum  control -surface  deflections: 

<5e,  deg .  10,  -25 

6,  deg .  ±15 

<3 

deg 


±30 
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TABLE  II  Aerodynamic  Characteristics,  Aircraft  F 


Cn 


!  T 


*».  «to*‘ 

-40 

-30 

-20 

-10 

10  , 

20 

30 

40 

-40 

-30 

-20 

-10 

0  10 

20 

30 

40 

-10 

0  547 

0  436 

0.324 

0  145 

0.000 

-0.1 S6  1 

-0  .120 

-0.431 

-0.542 

-0  04  2 

-0.042 

-0.042 

-0  020 

0  000  0  016 

0.036 

0  036 

0  036 

•5 

0.585 

0  464 

0  343 

0  158 

•0.143  j 

-0.320 

-0.441 

-0.561 

-0.042 

-0.043 

-0  044  -0  021 

|  0.016 

0.030 

0  037 

0  036 

0 

0  583 

0.467 

0  351 

0  167 

|  -0.126 

-0.309 

-0  425 

-0  540 

-0  039 

-0.041 

-0042 

-0  022 

0.015 

0  038 

0  037 

0  03s 

5 

0  572 

0.459 

0  346 

0.173 

:°:m6  j 

-0  284 

-0  397 

0.510 

017 

-0  038 

-0  039 

-0  022 

0  015 

0.027 

0  036 

0  "35 

10 

0  543 

0  443 

0  342 

0  180 

!  -0.102  1 

-0  26f> 

-0  367 

-0  401 

-0  035 

-0  034 

-0  034 

-0  021 

0  016 

0  032 

0  032 

0  012 

15 

0  506 

0.430 

0  353  i 

0.186 

!  .0.017  j 

-0  256  t 

-0  339 

-0  421 

-0  029 

-O  037 

-0  024 

-0  016 

0  012 

0  024  , 

0  026 

0  0?i» 

20 

'  0  471 

0  402 

0  .132 

0.109 

1  .0  004  ; 

0  221 

0  297 

0  370 

-0  005 

-0  003 

-0  001 

-0  008 

0  005 

0  008 

0  008 

0  *11  1 

25 

0  484  i 

0  390 

0  295  , 

0  161 

-0  051  J 

-0  187 

-0.282 

-0  376 

0  031 

0  028 

0  024 

0.007 

-0  006 

•0.015 

-0  010 

-0  004 

30 

0  544  j 

0.420 

0  295 

0  139 

-0.017  j 

•0  19b  . 

-0.323 

-0.447 

0  029 

0  032 

0  035 

0  024 

-0.021 

-0  029 

-0  020 

-0  01 1 

35 

0.618 

0  484 

0  350 

0.172 

-0.042  j 

-0  242 

-0  377 

0  511 

0.016 

0  025 

0  034 

0  028 

-0  027 

-0.033  ; 

-0  025 

-0  OK. 

40 

0  672  ! 

0.541 

0  410 

0  219 

•  0.097  , 

-0.285 

-0.410 

-0  547 

0.035 

0  032 

0.030 

0.025 

-0.021 

-0  031  , 

-0  033 

-0  036 

- . — 

» 

—  -  —  + 

...  i 

7 

—  . 

- - 

.  —  . 

45 

0  652  | 

0  532 

0  412 

0  254 

-0.126 

-0  .101 

-0.422 

-0  542 

0.04  0 

0  033 

0  020 

0  017 

-0  013 

-0.021 

■0  034 

-0  04 1'. 

50 

0.612  ; 

0  503 

0  394 

0  239 

!  .0.127  : 

-0.30C 

•0  416 

-0.525 

0  038 

0  016 

-0  005 

-0  015  : 

0  003 

0  004 

-0  017 

-0  036 

..  55 

0  618 

0  496 

0  373 

0.190 

-0.1197  ; 

-0  287 

-0.410 

-0  532 

0.034 

0.007 

-0.020  -0.039 

0  037 

0  020  . 

-0  007 

-0  031 

60 

0  670  . 

o.sn 

0  363 

0  153 

-0  051 

-0  261 

0  424 

-0.501 

0  044 

0.012 

-0  020 

-0  048 

0.040  0.013  j 

-0.018 

-0  050 

« 

1 

0.715  . 

0.574 

0.433 

0.120 

-0.049  [ 

-0.303 

0  444 

-0  585 

0  071 

0  042 

0  013 

-0.040  ; 

0  038 

-0  003  . 

-0  032 

-0  0*1 

70 

0.681  • 

0.508 

0  494 

0.174  ■ 

-0.075  1 

•  0  359 

-0  453 

-0  546 

0  072 

0  053 

0  033 

-0.023 

0  021 

-0  021  . 

-0  040 

-0  060 

75 

0.638 

0.505 

0  492 

0.247 

.0.107 

-0  376 

-0.451 

-0.524 

0  054 

0.03t 

0.022 

-0  008 

0.000 

-0.022 

-0  038 

-0  054 

80 

0.819  j 

0.551 

0  482 

0.289 

-0.143  ; 

-0  378 

-0.447 

-0.515  ' 

0  039 

0  022 

0  005 

-0.012  ! 

-0  007  |  -0  013  , 

-0.030 

-0.047 

L  '» 

0.613  1 

0.546 

0.476 

0.310 

-0.1C7  | 

-0  370  i 

-0.440 

-0.509 

0  037 

0.016 

-0.004 

0  014  j 

. 

-0.004 

-0.000  ! 

-0  023 

-0  040 

90 

0.620  j 

0.546 

0  471 

0  317  ! 

1  -0.184  ■ 

-0  355 1 

-0  430 

-0  504 

0.027 

0  009  |  -0  009  (  -0.015  ( 

j  0.008  |  0  002  | 

-0.017 

-0  011 

\ji.  deg 

0 

c 

m 

• 

or.  deg  \ 

I 

-40  | 

.... 

-30 

I 

-20  1 

-10 

0 

l 

10 

i 

20 

30 

40 

-40 

-30 

F 

-20 

-10 

10  ^  20  | 

30 

40 

-10 

-0.030  j 

-0.017 

-0.004  , 

0.001 

0  000 

-0  004 

0  004  ; 

0  017 

0  030 

0  119 

°.55 

0  181 

0.209 

0  32*1  0.291 

0  240 

0  184 

0  1 t92 

-5 

-0.011  I 

-0  003 

0.006 

0.003 

-0.005 

-0.007 

-0.004 

0  000 

0  330 

0  240 

0  155 

0  189 

0 

173  0170 

0.186  , 

0  261 

0  3  i». 

0 

0.024  , 

0.023 

0  022 

0.000 

-0.012 

-0.022 

-0  024 

-0.025 

0  420 

0.200 

0  113  1  0  081 

0.003  0  063  0  117  j 

0  269 

0  420 

5 

0.051 

0.046 

0  040  ' 

0.017  ] 

-0.02! 

-0.039 

-0.045 

-0.050 

0  418 

0  225 

0.032 

-0  031 

-0.037  -0  050 

0.026 

0  222 

0  41h 

10 

0.076 

0.060 

0  044 

0.020  , 

-0.023 

-0.046 

-0  062 

-0 .078 

0  418 

0.160 

-0  090 

-0  142 

-0.148  -0  182 

-0.090 

0.164 

0  41b 

15 

0.096 

0  068 

0  040 

0.020 

-0.024 

-0  043 

-0.071 

-0.099 

0  378 

0  074  ,  -0.229 

-0.222 

-0  218  -0  220 

-0.215  0.081 

0  37K 

20 

0.096 

0.063 

0  033 

0  019  , 

,  -9.023 

-0  018 

-0.070 

-0.101 

0  326 

0  013 

-0  299 

-0  307 

-0.284  -0  309 

-0.311 

0.003 

0  326 

i* 

0.075 

0.047 

0  019 

0  OH.  ■ 

-9  019 

-0  031 

-0  059 

-0.087 

0  386 

0  030 

-0  320 

-0.259 

-0.401  -0  301 

-0.301 

0  013 

0  386 

!  30 

o.ors 

0.039 

0  012 

0.007  : 

1  -0.010 

-0.015 

-0.042 

-0.068 

0  491 

0  083 

-0  324 

-0.343 

-0.531  j  -0.437 

-0  347 

0.072 

0  491 

f  35 

0.051 

0.033 

0.015 

-0.00*- 

0.000 

-0  012 

-0.030 

-0.048 

0  535 

0  124  |  -0.287 

-0.377 

-0.579  j  -0  412  -0,337  ; 

0.099 

0.535 

t  40 

0.038 

0-023 

0  008 

-0.007  . 

0.001 

-0  017 

-0  032 

-0.017 

0  48t 

0  121 

-0  244 

-0.493 

-0.603  -0.439 

-0.340  , 

0.073 

0  466 

45 

0.060 

0.037 

0  014 

-0.005  1 

-0.006 

-0  023 

-0.046 

-0.069 

0.247 

0.037 

-0  174 

-0  115 

-0.617  j  -0.550 

-0.265 

-0.009 

0  247 

50 

0  091 

0.055 

0  038 

0.002  ‘ 

;  .0  011 

-0  030 

-0  057 

-0  083 

0  073 

-0.049 

-0  170 

-0  444 1 

-0  026  -0  507 

-0.091 

-0  000 

0  on 

S3 

0  089 

0.070 

0  050 

0.0 1 8  ! 

-0.019 

-0  045  , 

-0.005 

-0  084  ’ 

0  024 

.0  in 

-0.249 

.0  515  I 

-0  647  |  .0  410 

-0  102 

-0  079 

C  024 

00 

0.087 

0.070 

0  052 

0.076 

-0.026 

-0  051 

-0.009 

-0  088  ' 

-0  105 

-0  299 

-0  412 

-0.631 

-0  703  -0  402 

-0.333  , 

-0  259 

-1.185 

05 

0.091 

0.072 

0  052 

n  075 

|  -0  028 

-0  052 

-0  070 

-0  088 

-0  494 

-0.550 

-0  871 

-0.705 

•  0  605  !  -0  648 

0  470 

-0  hi 

-3  494 

70 

0  092 

0.072 

0  052 

0  027  1  1 

-0.027 

-0  053 

-0  072 

-0.030 

-9  719 

-0.773 

-0.026 

-0  068 

-0.953  \  -0.012  j  -0  706 

-0  713 

-0.719 

75 

0.093 

0  073 

0  052 

0.026 

-0  016 

-0  052 

-0  071 

•0  033 

-0  868 

-0.975 

-1.082 

-1  002 

-1 

130  J  0  990 

-0.999  j 

-0.(34  . 

-0  86H 

00 

o  o«5 ; 

o  on 

0  048 

0.027 

0.000 

0  050 

-0.071 

•  n  078 

-1  000 

-1.100 

-1.271 

-1  294  | 

-1.328  :  -1  201 

-1  SOS 

•  ^ 

-1.154 

-t  000 

05 

0.085  I 

0  072 

0  049t 

0.029 

! 

-0  006 

0.061 

-0  074 

-0-0i7 

-1  135 

-1.351 

-1  567 

-1  541  1 

-1 

619  k  -1.492 

-1  329 

-1  135 

90 

O.DH  j 

0.073 

0  050 

0.032 

* 

-0.027 

-0  052 

-0.073 

-0  098  ’ 

-1  274 

-1  407 

-1.700 

-i.*ii  | 

-1.974  j  -1  043 

,6i2! 

-1  473  | 

-l  274 
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TABLE  II,  concluded 


.... 

nr. 

deg 

ex 

IT 

Cv  , 

5> 

Per  deg 

7,7 

\ 

per  deg 

C"V 

pci  deg 

7 

S 

per  rad 

7~ 

V 

per  deg 

s  • 

r 

per  dec 

-0.0090 

0.8500 

0  0001 

-0.0121 

-0.0264 

-26  040 

0.0035 

7  “7 

r5 

-0.0250 

0.4600 

0  0035 

-0.0150 

-0.0297 

•  26.040 

0.0034 

-0  0013 

0 

-0  0300 

0.0600 

0.0050 

-0.0150 

-0.0303 

-20.040 

0.0032 

-0  0013 

5 

-0.0286 

-0.3200 

0.0043 

-0.0140 

-0.0300 

-14.420 

0.003li-0.00l2 

10 

0  0000 

-0.7300 

0.0035 

-0.0143  -0.0305 

-22.790 

0.0029 

-0.0013 

15 

0  0451 

-1.1300 

0.0026 

-0.0154 

-0.0312 

-26.230 

0.0030 

-0.0013 

20 

0.0131 

-1.5300 

0.0017 

-0.01*4 

-0.320 

-29  670 

0.0032 

-0.0013 

0.0991 

-1.9200 

0.0006 

-0.0207 

-0.0344 

-33.700 

0.0033 

-0.0013 

j? 

0.0920 

-2. 3300| -0.0003 

-0.0241 

-0.0370 

-37.720 

0.0032 

-0.0013 

t** 

0.0708 

-2  6500 

-0.0012 

-0  0250 

-0  0361 

-42.900 

0  0029 

-0.0012 

40 

0.0465 

-1  7470 

-0.0011 

-0.0220 

-0.0316 

-44  690 

0.0025 

-0.0009 

0.0397 

-1.6*96 

-0  0035 

-0.0169 

-0  0247 

-41.810 

0  0019 

-0  0007 

1  5° 

0.0334 

-1  7054 

-0  0043 

-0.0144 

-0.0174 

-19000 

0.0019 

-0.0006 

!  55 

0.0392 

-1.7492 

-0  0048 

-0.0120 

-0.0114 

2  000 

0.0041 

-0.0006 

I  60 

0  0391 

-1.7680 

-0.0045 

-0.0103 

-0  0011 

-7  000 

0.0016 

; 

I  65 

0.0344 

-1.8142 

-0.0047 

-0.0093 

-0  0034 

-30.000 

-0.0005 

|  TO  1  0.0350 

•  1.9020 

-0.0032 

-0.0105 

-0  0040 

-11.000 

0.0004 

J 

1  75 

0  0395 

-1.9490 

-0.0053 

•0.0106 

-0.0040 

-4.000 

0.0005 

! 

80 

0.0407 

-1.9634 

-0.0052 

-0  0082 

-0.0040 

-3.000 

-0.0000 

:: 

|  85]  0.0412 

•1.9690 

-0.0050 

-0.0071 

-0  0040 

-17000 

-0.0005 

L90 

0.M12 

-1.9690] -0.0061 

-0.0077 

-0.0040 

-24  000 

-0.0008 

J 

•Kis  ’ 

‘K  per  deg  per  d»  * 
12 1  0  OOOS^O  0002 
12 1  0  0005  -0  0002 
•2^  0  0011. -0  0002 
>2 '  0  0006  -0  0001 
>2*  0  0006  -0  0001 
0.0001,  0  0010  0  0000 
0  0001 
,  0.0001 
0  0002 

0  0003 

- - 

0  0003 


S' 


per  deg  per  rad  per  rad  per  rad  per  rad  per  'ad'pcr  rad 


ft 


—4 


0004,  0  0007 

0003 j  0.0002 

0  0001 1-0.0013 
-f-  -4 

0  0001^0  0035 

0  OOOOj-O  0035 

0  0000 i-0  0027 

•0.0025, 

0  0024 

.-•0  0024 


0  0008]  0  030 

-0  010  ' 

0  140 

0  120  :  -0  170  0  040 

0  0008  0  060 

-0  010 

O'  190 

&  130  -0  170^  0  040 

0  160  !  -0  170  |  0.040 

0  0008 ,  0  120 

-0  010 

0  190  . 

0  OOOg|  0  190 

-0  010 

0  160 

0  1  30  -  0  170  .  0  090 

■  --  .  -  . 

- *  1 

0.0009'  0  230 
0.0008^  0.240 

0  00081  0  230 

-0  010 

0  180 

0.010  ;  -0  180^  0  130 

0  000  ,  -0  2  20  0.220 

0.000 

0  180 

0  010 

0.160 

0  430  -0  280  ]  0  310 

0  0009]  0  260 

0  190 

0  180  , 

1  050  -0  270  0  480 

0  0009  j  0  290 

0  360  : 

0  200 

1  200  -0  280  0  640 

0  0009 |  0  290 

0  580 

0  380 , 

0  790  -  0  280  1  160 

0  00091  0  500 

0  400 

0  550 

0  230  j  -0.1*0  ,  1  590 

0  0009 i  1  230 

0.260  * 

0  800: 

0  180  -0  080  0  990 

0  0007  1  700  0  190  -0  570  0  780  0  160  0  350 

0.00061  1.540  0  140  -0  450  2  610  f  0  990  ?  0  250 

-0  310  !  -<>  270  2  270  '  0  820  ;  0  130 


O.OOO'yO  140 
0  0005  -1180  -0  470  i-0  ISO!  0.490T  0.000 
0  0004  ;  -0  090  -  0  050  -0  100  ,  -0  180  -  0  1 10 


0 .000 5 j  0  640  ; 
0  00061  0  580  1 

4  - 

0  0005]  0  610  ; 
-0  0002!  0  730  ' 


-0  150  I  -0  ioo-  -0  010  ;  -0  no 

0  090  I  -O  HO 


0  040  I  -0  133 
-0  040  j  -0  140 
-0  050  |  -0.150 


0  090  ^ 

0  100  i  -0  110 


1  j  •»  J 1' 
)[-0  111 


0  030  ; 
0  000  1 

- - H 

0.010 

1 


0  020 
0.010 
0.010 


1 
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The  general  form  of  a  bi-cubic  spline  is 
4  4 

S(a,B)  =  l  l  C..[a-A(st  )]1"1[6- B(ft R)]3'1  (3.2.52) 

i=l  j=l  a  y 

where  A( 2,^)  and  B(£g)  are  the  values  for  a.  and  8,  respectively,  at  the 
lower  left  corner  of  the  rectangle  of  values  for  a  and  8  which  contains 
the  input  set  (a, 8).  S  is  the  value  of  the  particular  coefficient  at  the 
input  set  (a, 8).  Library  routines  (IMSL)  are  used  to  generate  the  C^, 
and  to  compute  S  and  its  parti als. 

Use  of  splines  may  appear  to  be  introducing  an  overly-sophisticated 
approach  for  a  data  base  as  relatively  simple  as  the  aircraft  F  model,  but 
when  the  more  complex  data  bases  such  as  that  for  the  F-4  are  added,  the 
value  of  the  spline  approximations  should  be  more  appreciated.  To  assure 
a  good  spline  fit,  the  location  of  the  junction  points,  or  knots,  is  an 
important  factor.  When  this  and  similar  considerations  are  efficiently 
dealt  with,  the  method  of  splines  becomes  an  efficient  modeling  tool.  In 
our  applications,  the  knots  are  placed  at  each  data  point.  More  back¬ 
ground  on  the  theory  and  uses  of  splines  may  be  found  in  Ahlberg,  Nil  son 
and  Walsh  (1967). 

Splines  by  their  nature  allow  for  accurate  modeling  of  the  partial 
derivatives  of  the  dependent  variable  (aerodynamic  coefficients)  with 
respect  to  the  independent  variable(s)  (a  and/or  e).  This  fact  is  impor¬ 
tant  in  the  application  here,  because  these  partials  are  required  in 
jrder  to  evaluate  the  equilibrium  and  bifurcation  surfaces;  and  the  accuracy 
*  •"('ir  values  may  often  be  critical  to  efficiently  starting  up  the 
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numerical  solution.  Whereas  curves  which  result  from  a  less  rigorous 
patching -tog ether  of  polynomials  may  well  experience  severe  discontinuity 
problems  for  their  derivatives  at  the  boundary  points,  spline-produced 
curves  by  their  nature  (i.e.,  the  constraint  that  the  second  derivatives 
be  equal  at  each  knot  which  "shares"  two  splines)  have  no  such  problems. 

The  splines  used  to  model  aircraft  F‘ s  aerodynamic  data  are  "natural" 
cubic  (or  bi-cubic)  splines*,  and  the  knots  are  specified  to  be  located 
at  each  data  point.  Of  the  22  coefficients  which  are  used  by  aircraft  F, 
all  but  four  are  functions  of  a  and  are  thus  modeled  by  the  one-dimen¬ 
sional  cubic  splines.  C^,  Cn,  C^,  and  only  are  functions  both  of  a  and 
P  and  are  modeled  by  bi-cubic  splines.  Cubic  spline  plots  are  shown  in 
Figs.  3.1  to  3.6.  These  are  plots  of  each  of  the  18  a-dependent  coeffi¬ 
cients,  and  their  respective  partial s  with  respect  to  a  (in  degrees). 

The  plot  of  the  basic  function  is  the  line  containing  the  x's,  the  latter 
representing  the  data  points  (hence,  also  the  knots)  obtained  from  Moore, 
Anglin  (1971).  Both  the  coefficients  and  their  partial s  with  respect  to 
a  are  computed  and  plotted  in  dimensionless  units.  The  units  for  a  (abscissa) 
are  degrees.  The  derivative  with  respect  to  a  is  plotted  as  a  clean,  solid 
line.  Notice  how  smooth  the  derivative  curves  are.  The  large  change  in 
shape  of  most  of  the  curves  beginning  in  the  range  a = 35°  to  50°  (.69  to  .91 
radians)  is  due  principally  to  the  loss  of  rudder  control  effectiveness  for 
a  >50°.  See,  in  particular,  the  plots  of  the  lateral  mode  coefficients 
appearing  Jn  Figs.  3.3  to  3.6. 

The  bi-cubic  spline,  two-dimensional,  plots  of  a  representative 

function  of  both  a  and  e,  Cm(a,8),  are  presented  in  Figs.  3.7  and  3.8. 

*A  natural  cubic  spline  is  one  whose  second  derivatives  are  zero  at  the 
end  points. 
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Fig.  3.7  shows  Cm  and  its  B-deri vati ve,  Cm  ,  both  of  which  are  in  dimen¬ 
sionless  units,  plotted  against  sideslip  angle  e,  in  degrees.  The  two 
plots  shown  in  Fig.  3.7  are  for  two  different  values  of  angl e-of-attack 
a,  65°  and  32.5°  respectively.  The  latter  value  for  a  represents  an  in¬ 
terpolation  on  the  data.  By  rotating  the  projection  plane  depicted  in 
Fig.  3.7  by  90°,  one  obtains  plots  of  C  (a,e)  and  its  a-derivative,  C  , 

a 

versus  a,  shown  in  Fig.  3.8  .  The  units  are  as  in  Fig.  3.7  .  The  two 
values  of  B  for  which  the  Fig.  3.8  plots  are  made,  +25°,  are  the  interpola¬ 
tion  values  for  6. 

A  final  note  on  the  use  of  spline  function  approximations  by  BACTM. 
Because  the  spline  package  had  already  been  implemented  to  represent  the 
aero  data,  it  was  also  decided  to  use  splines,  and  the  same  routines, 
in  evaluating  numerical  derivative,  e.g. ,  to  generate  elements  of 
[8f/3x]  which  are  difficult  to  obtain  analytically.  This  duality  of 
function  adds  to  the  efficiency  of  BACTM.  See  Section  2.1.2  for  a  dis¬ 
cussion  of  the  numerical  differentiation  algorithm. 
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3.3  Entry  into  Spin;  Explanation  of  Equilibrium  Surface  Plots 

Experience  shows  that  a  very  complete  set  of  aerodynamics  is  required 
to  represent  properly  the  extremely  complex  aircraft  motions  which  arise 
during  stall/spin  flight  conditions.  In  this  section,  we  shall  explore 
the  spin  entry  behavior  of  aircraft  F.  This  model  does  not  possess  all  of 
the  desirable  features  for  undertaking  a  thorough,  realistic  study  of  spin 
behavior.  However,  it  has  been  an  excellent  model  in  providing  a  small 
yet.  adequate  basis  about  which  to  construct  the  BACTM  Spin  Analysis  system; 
and  for  providing  much  insight  into  the  nature  of  developed  spin  motion  in 
particular. 

The  results  to  be  shown  in  this  section,  however,  will  reveal  that 
careful  study  of  the  equilibrium  surfaces  generated  by  BACTM  is  required 
in  order  to  understand  the  subtleties  of  nonlinear,  high-a  aircraft 
motion.  Interpretation  of  equilibrium  surfaces  is  a  case  in  point.  It 
can  be  rigorously  demonstrated  that  motions  originating  in  the  vicinity 
of  a  "stable"  branch  (such  branches  are  labeled  with  an  "S"  in  the  figures) 
will  ultimately  arrive  at  a  point  on  that  branch,  thereby  achieving  a 
condition  of  dynamic  equilibrium  (i.e.,  x  =  0).  See  Section  3.2  and  Chapter  II 
for  further  discussion.  Chapter  II  contains  a  bibliography  referencing  detailed 
proofs.  The  only  problem,  then,  for  regions  dominated  by  a  stable  branch 
is  determining  quantitatively  the  neighborhood  from  which  all  motions  lead 
to  that  branch  (this  neighborhood  is  called  a  "domain  of  attraction"). 

Very  often,  a  boundary  to  the  domain  of  attraction  is  a  "simple"  un¬ 
stable  branch,  e.g. ,  one  in  which  there  is  one  positive  real  root--this 
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branch  is  designated  by  a  "U".  In  Fig.  3.9  we  can  see  an  example  of 
this  situation.  The  point  A  in  this  figure  is  centered  on  a  stable  branch 
in  the  equilibrium  flat  spin  region.  This  branch  is  bordered  above  and 
below  by  a  U-branch.  Thus,  motions  initiated  with  a  contrul  setting  close 
to  that  of  point  A,  and  with  initial  conditions  inside  the  two  U-branches, 
will  return  to  the  S-branch.  Units  for  this  and  subsequent  figures  are  degrees. 

Before  proceeding  further,  we  will  now  explain  some  of  the  details 
of  Fig.  3.9,  which  is  a  typical  plot  of  an  equilibrium  surface.  The 
branches  shown  represent  loci  of  points  for  which  f(x,5)  =  0,  generated  by 
holding  two  of  the  elements  of  6  at  given  values,  and  varying  the  third 
(see  Chapter  II).  The  letters  on  these  branches  provide  local  stability 

X.  L. 

information,  as  follows:  for  the  ntn  order  system  x  =  f(x,6),  the  local 
stability  information  is  obtained  from  the  eigenvalues  of  the  square 
matrix  [9f/9x],  which  is  part  of  the  first  order  (linear)  term  in  the 
polynomial  expansion  of  f.  This  matrix  is  of  size  nxn,  and  always  yields 
n  eigenvalues,  which  may  either  be  real  or  in  complex  pairs  (since 
[9f/9x]  has  real  elements).  If  all  n  of  the  eigenvalues  have  negative 
real  parts,  the  equilibrium  point  is  a  stable  one,  designated  by  S.  Any 
other  situation  results  in  an  equilibrium  point  which  is  locally  unstable. 
Several  unstable  cases  are  now  outlined:  if  there  are  (n-1)  eigenvalues 
with  negative  real  parts,  and  the  remaining  one  is  positive  (this  one  is 
necessarily  a  real  eigenvalue  as  complex  eigenvalues  come  in  pairs),  the 
point  is  designated  by  a  U;  if  (n-2)  have  negative  real  parts  and  a  com¬ 
plex  pair  has  a  positive  real  part,  the  point  is  designated  by  L.  How- 

* 

ever,  if  there  are  (n-2)  "stable"  eigenvalues  and  two  real,  positive  ones, 
then  the  point  has  two  simple  U's,  or  is  designated  by  the  symbol  A,  as 
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seen  in  the  table  at  the  top  of  Fig.  3.9.  Except  for  the  one,  unique 
stable  case  (S),  the  table  gives  symbolic  correspondences  for  all  com¬ 
binations  of  unstable  eigenvalues;  that  is,  eigenvalues  with  positive 
real  parts. 

To  summarize  to  this  point,  when  U-branches  have  been  found  near 
S-branches,  they  typically  outline  much  of  the  domain  of  attraction  to 
that  stable  branch. 

As  for  all  other  (unstable)  branches,  there  are  two  major  points  of 
interest  with  regards  to  aircraft  F  equilibrium  results: 

1)  The  stability  analysis  is  local,  and  based  only  on  first  order, 

or  linear,  information.  This  means  that  at  a  U-branch,  say  the  one  including 

point  E  in  Fig.  3.9,  the  motion  will  diverge  locally;  in  fact,  for  an^ 

unstable  point,  any  point  but  an  S-point,  the  motion  diverges  locally. 

But  we  can  say  more  regarding  the  point  E  case:  if  the  motion  starts  on 

the  S-branch  side  of  the  U-branch,  it  will  ultimately  come  to  equilibrium 

on  the  S-branch.  This  is  a  global  result  however,  and  is  obtained  not 

from  the  information  provided  by  [3f/3x]  at  point  E,  but  from  a  more  global 

knowledge  of  the  equilibrium  surface--that  is,  we  know  that  the  S-branch 

"attractor"  exists.  If  the  motion  starts  on  the  other  side*  of  the  U- 

branch,  it  is  unclear  from  the  equilibrium  surface  plot  what  will  ultimately 

happen.  The  divergence  feature  applies  only  in  a  local  sense  and 

guarantees  only  that  there  will  be  no  equilibrium  near  the  U-branch. 

Globally,  the  motion  is  attracted  elsewhere  from  the  U-branch,  either  to 

a  distant  S-branch  (not  shown),  or  to  a  limit  cycle,  perhaps  the  one 

governed  by  the  L-branch  above  point  E  in  Fig.  3.9a;  or  finally,  to  an 

♦Strictly  speaking,  motions  originating  exactly  on  the  U-branch  or  any 
other  unstable  branch,  will  remain  there;  however,  the  smallest  pertur¬ 
bation  will  induce  divergent  behavior,  which  is  what  happens  in  the  real 
world. 
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essentially  erratic  motion  not  characterized  as  either  a  stable  equi¬ 
librium  (x  =  0)  or  a  limit  cycle  (x(t  +  T)  =  x(T) ,  0  <  T  <  °°) .  Thus,  local 
instability  is  not  necessarily  global  instability. 

2)  L-branches  are  of  particular  interest,  because  they  are  very 
prevalent  in  high-a  flight  regimes,  and  particularly  so  in  the  spin 
entry  region  depicted  in  Fig.  3.9  (the  lower  branch  which  includes  points 
B,  C,  D,  E).  As  defined  above,  the  L-branches  define  unstable  equilibrium 
points  at  which  [3f/9x],  the  Jacobian  matrix,  has  (n-2)  "stable"  eigen¬ 
values  and  one  unstable  complex  pair.  Locally,  then,  this  indicates 
oscillatory  divergence.  Globally,  however,  it  is  quite  possible  that 
this  ivergence  is  in  actuality  growth  to  a  stable  limit  cycle  (hence, 
the  L-symbol  designation).  The  existence  and  stability  of  limi't  cycles 
is  not  explicitly  obtained  from  the  equilibrium  surfaces  (we  are  currently 
developing  an  algorithm  to  do  this  analysis),  but  stable  limit  cycles  in 
our  experience  are  always  associated  with  an  L-branch.  In  a  sense,  then, 
stable  limit  cycles  may  be  regarded  as  stable  equilibria,  unless  ampli¬ 
tude  variations  are  excessive.  At  any  rate,  there  are  domains  of  attrac¬ 
tion  to  stable  limit  cycles,  although  their  boundaries  are  not  as  easily 
computed  as  are  those  for  the  proper  equilibrium  points  (S-branches) . 

Also,  there  are  known  to  exist  unstable  limit  cycles.  In  analogy  to  the 
U-branches,  if  a  motion  could  begin  exactly  on  an  unstable  limit  cycle, 
and  be  free  of  random  disturbances,  it  would  remain  on  this  cycle.  Typically, 
of  course,  the  motion  readily  diverges,  usually  to  a  stable  limit  cycle. 
Unstable  limit  cycles  are  very  difficult  to  isolate,  and  there  exists  no 
known  algorithm  of  sufficient  generality  which  can  compute  their  location. 

They  do  exist  however,  and,  as  the  continuation  parameter  (6r  in  Fig.  3.9) 
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changes,  are  capable  of  annihilating  a  stable  limit  cycle.  This  phenomenon 
is  called  a  global  bifurcation,  and  may  explain  the  limit  point  below 
point  B  in  Fig.  3.9a.  Here,  as  one  decreases  6r  from  point  D,  the  equi¬ 
librium  curves  tell  us  to  expect  a  Hopf  bifurcation  (S-to-L  transition 
on  the  branch)  to  a  limit  cycle.  Such  limit  cycles  are  typically  stable, 
at  least  when  in  the  region  close  to  the  S  branch,  but  again,  this  in¬ 
formation  cannot  be  obtained  from  Fig.  3.9.  As  6r  decreases,  however, 
this  limit  cycle  is  annihilated;  we  cannot  indicate  for  sure  exactly 
where  (although  the  figure  indicates  approximately  -21°,  the  limit  cycle 
motion  is  not  a  local  motion',  and  the  limit  cycle  may  be  annihilated 
at  a  very  different  value,  if  at  all),  but,  given  that  this  occurs,  it 
is  quite  likely  that  the  motion  will  be  attracted  to  a  stable  limit  cycle 
governed  by  the  segment  of  the  equilibrium  branch  including  point  B.  This 
hypothetical  control  sequence  outlines  how  a  relatively  "clean"  rolling 
motion  at  point  D  (6a =  15°  in  this  example,  recall)  is  corrupted  into 
buffeting  and  oscillatory  motion  which  conceivably,  as  6r  is  decreased 
further,  undergoes  a  jump  to  oscillatory,  steep  spin  conditions  (point  B). 
From  analysis  of  time  history  results,  discussed  in  more  detail  later, 
there  is  a  strongly  stable--i.e. ,  possessing  a  large  domain  of  attraction-- 
limit  cycle  in  the  vicinity  of  point  B.  This  situation  precludes  transi¬ 
tion  from  steep  spin  to  the  flat  spin  equilibrium  at  point  A,  using  6r 
alone. 

Summarizing  unstable  branches,  then,  the  major  point  is  that  in  a 
local  sense  all  motions  diverge  from  them;  however,  limit  cycles  typically 
exist  about  a  subset  of  these,  the  L-branches,  and  these  themselves  may 
be  stable  (attractors)  or  unstable.  When  a  stable  limit  cycle  exists. 
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it  is  usually  centered  on  the  L-branch;  however,  the  motion  will  never 
decay  to  the  L-branch  itself.  In  spite  of  this,  it  is  often  convenient 
to  consider  the  stable  limit  cycle  as  an  equilibrium  condition,  which  may 
be  disrupted  as  a  neighboring  unstable  limit  cycle  converges  upon  it, 
under  the  influence  of  changing  control  parameters. 

A  final  note  on  limit  cycles:  they  are  a  distinctly  nonlinear  phe¬ 
nomenon,  lacking  many  properties  of  the  linear  oscillator.  In  the  latter, 
an  incremental  change  in  initial  conditions  results  in  an  incrementally 
different,  yet  stable,  orbit.  Such  a  change  applied  in  a  limit  cycle 
region  will  only  produce  a  temporary  perturbation,  followed  by  decay  to 
the  original  limit  cycle.  That  is,  in  a  given  subspace,  there  are  only 
a  countable  number  of  possible  limit  cycles,  while  a  continuum  of  ampli¬ 
tudes  is  possible  for  the  linear  oscillator. 

Returning  to  the  aircraft  F  spin  entry  case,  Fig.  3.9  is  a  repre¬ 
sentative  situation  (the  insensitivity  to  6e  of  the  equilibrium  surfaces 
in  this  flight  regime  will  be  indicated  later;  thus,  6e  =  0°  is  selected), 
in  that  the  flat,  equilibrium  spin  region,  the  upper  branch  in  Fig.  3.9a, 
has  a  relatively  small  stable  equilibrium  region.  Further,  the  steep  and 
intermediate  spin  regions,  the  upper  part  of  the  lower  branch,  as  well  as 
most  of  the  other  non-equilibrium  regions,  are  characterized  by  oscilla¬ 
tory  behavior,  with  limit  cycle  motion  quite  likely.  As  stated  above, 
several  of  the  limit  cycle  regions  exert  strong  attraction  on  neighboring 
motions,  making  transition  to  flat  spin  very  difficult.  At  any  rate,  our 
experience  with  aircraft  F  indicates  that  oscillatory  motions  at  inter¬ 
mediate  angles  of  attack  (30°  s  a  s  65°)  is  a  general  feature  for  any  control 
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setting;  this  observation  consequently  adds  to  our  desire  to  study  in 
detail  the  behavior  and  stability  of  limit  cycle  motion  in  general,  and 
the  stall/post-stall/spin  entry  flight  regime  in  particular,  because 
of  its  oscillatory  character.  A  more  useful  study  of  spin  entry  would 
involve  using  the  F-4  model,  which  is  more  realistic  and  representative 
in  this  flight  regime. 

In  addition  to  the  above  situation,  Bihrle  (1976)  has  noted  (and  it 
is  verified  here)  that  ensuing  high-a  motion  is  extremely  sensitive  both 
to  control  sequencing  and  to  relatively  small  variations  in  the  initial 
conditions.  Results  presented  here  provide  the  basis  for  the  above 
observations  and  are  typical  of  systems  with  bifurcations. 

In  the  following  discussion,  we  have  concentrated  on  right  pro-spin 
motions;  that  is,  yaw  rate  (r)  is  positive,  due  to  negative  rudder  (fir) 
deflection,  and  aileron  (fia)  is  positive.  For  convenience,  then,  we  shall 
def i ne 

-SPIN  ~  Ua, fie, fir)  =  (15,-21,-25)  degrees 

Because  of  nonlinear  dynamics,  this  set  of  controls  corresponds  to  more  than 
one  equilibrium  point.  However,  associated  with  fi<jpjN  is  a  stable,  de¬ 
veloped  flat  spin  equilibrium  state  for  aircraft  F, 

-SPIN  —  ^  P » ^  » 6 ,  V ,  9 ,  <j> ) 

=  (30, -4., 100., 73. 5, -3. ,443,-16.6,-2.29) 

The  units  are  degrees  and  feet  per  second.  These  quantities  will  be  used 
for  reference  in  the  following  sections. 

Finally,  it  should  be  noted  that  many  of  the  equilibrium  surfaces 
presented  in  this  chapter  may  not  contain  all  of  the  possible  equilibrium 
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branches  for  the  given  control  region.  This  is  not  necessarily  a  minor 
detail,  because  some  of  the  branches  not  shown  may  possibly  represent 
stable  attractors  for  certain  limit  cycle  motions. 

3.3.1  High-a  Motions  to  Stall,  Aircraft  F 

Fig.  3.10  shows  equilibrium  surfaces  in  the  (r,a,p)-6e  plane  for 
aircraft  F,  centered  on  the  trim  state  with  neutral  controls.  As  will  be 
the  case  unless  otherwise  specified,  velocity  (V)  is  fixed  at  600  fps,  and 
gravity  is  assumed  negligible  (i.e.,  the  non-spin  set  of  dynamic  equations 
is  used).  Note  that  there  are  regions  (5e>0)  where  five  equilibrium 
solutions  exist,  and  three  of  these  are  stable.  As  6e  goes  negative, 
pitch-up  occurs,  signified  by  the  growing  values  in  a.  From  6e  = -15°  the 
stable  branch  changes  into  a  limit  cycle  branch;  this  signifies  the  onset 
of  wing  rock  behavior  and  pre-stall  buffeting.  Note  in  particular  that 
on  the  roll-rate  plot,  there  is  a  region  in  6e  > 0  for  which  no  solutions 
are  shown.  At  the  boundary  points  to  this  region,  a  has  reached  its 
minimum  value  of  -10°,  and  no  aerodata  were  available  below  this  value. 

Fig.  3.11a  shows  a  time  history  plot  in  which  V  is  free  to  vary  and  6e  is 
increased  from  0°  to  -9°,  then  -17°,  then  -20°.  Note  that  when  6e  is 
-17°  and  -20°,  there  is  evidence  of  large  longitudinal  oscillations,  but 
little  wing  rock  motion.  Fig.  3.11b  shows  a  different  5e  time  sequence 
than  3.11a,  with  velocity  allowed  to  vary,  and  basically  similar  results. 

If  velocity  is  kept  constant,  buffeting  and  wing  rock  activity  are  very 
prominent  (Fig.  3.11c)  as  will  be  shown  in  Sec.  3.3.3,  where  Figs.  3.10 
and  3.11  are  discussed  further. 

We  consider  now  a  maneuver  in  which  both  heading  change  and  high-g 
pullup  is  accomplished  by  pulling  the  stick  back  sharply  (high  negative  5e), 


and  then  over  (high  aileron).  The  5e  maneuver  can  be  traced  in  Figure  3.10, 
and  Figure  3.12  shows  the  6a  direction  from  Point  A  in  Figure  3.10.  It  is 
seen  in  Figure  3.12  that,  for  6e  =  -11°,  beyond  6a  =  15°  there  is  a  jump, 
or  a  Hopf  Bifurcation,  to  a  limit  cycle  with  high  values  of  (r,i,p). 

The  time  history  presented  in  Figure  3.13  confirms  these  conclusions, 
although  the  transition  to  post-stall  divergence  occurring  at  30  sec.  is 
aided  by  introducing  negative  Sr  as  well.  Note  that  the  presence  of  nonzero 
6r  causes  the  initiation  of  spin-like  behavior  in  the  post-stall  motions. 

The  DSPIN  parameter,  a  normalized  metric  for  ||(x,6)  -  (*spin>§spin)  II 2  > 
shows  a  return  to  spin  conditions.  The  oscillatory  nature  of  the  developing 
spin  is  typical  of  aircraft  F  behavior  in  transitioning  from  trim  to  hi gh-cx 
flight  regimes. 

The  equilibrium  surfaces  in  Fig.  3.14  show  the  rudder  effect  on  the 
type  of  equilibrium  condition  which  produces  aircraft  behavior  similar  to 
that  of  the  maneuver  discussed  above.  In  this  instance,  6e =  -20°  and 
6a  =  0°,  while  6r  varies.  The  segment  of  the  L-branch  (limit  cycle  beha¬ 
vior  expected)  near  (p,r)  =  0,  6r=0°,  actually  exhibits  very  mild  unstable 
growth  to  oscillatory  behavior  (very  large  characteristic  or  response  time) 
when  the  corresponding  time  hi  story  run  is  made  (figure  not  shown);  but, 
as  6r  moves  from  0°  to  negative  values,  the  response  quickens  greatly,  and 
the  oscillations  grow.  A  similar  effect  is  noted  in  the  trajectories  of 
Fig.  3.13. 

3.3.2  Non-Spin  Equilibrium  Surfaces  (Aircraft  F) 

As  above,  these  results  were  generated  from  the  non-spin  set  of 


dynamic  equations,  which  assume  zero  gravity  and  constant  velocity  (V). 
Except  where  stated  otherwise,  V  is  assumed  to  be  600  fps.  It  turns  out. 
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at  least  for  aircraft  F,  that  the  non-spin  set  can  also  be  used  in  the 
SDin  regime  without  destroying  essential  features.  This  represents  a 
significant  result  for  first-cut  analysis;  the  numerical  results,  however, 
do  differ,  since  gravity  does  play  a  role  in  spin  behavior.  However, 
gravity  has  negligible  influence  on  equilibrium  behavior  in  non-spin, 
high-a  regimes.  Moreover,  steady  state  values  for  e  and  tp  do  not  exist 
in  the  high-a  regimes.  Hacker  and  Oprisiu  (1974)  show  that  the  effects 
of  gravity  in  roll  coupling  may  be  taken  into  account  by  a  perturbation 
analysis. 

The  projections  of  the  equilibrium  manifold  on  different  planes  in 
the  state  and  control  space  are  shown  in  Figures  3.10,  3.12,  3.14,  3.15 
and  3.16.  In  each  of  these  runs,  the  control  is  extended  from  x  =  0  over 
as  much  of  its  allowable  range  as  possible,  with  the  remaining  two  controls 
fixed  at  zero.  The  elevator  case  has  been  discussed  in  Section  3.3.1. 

As  will  be  seen,  the  rudder  introduces  the  most  dramatic  changes  in  equilibrium 
conditions  in  this  region  (Figures  3.14  and  3.15).  In  fact,  the  aileron 
behaves  as  an  almost  purely  linear  control  over  its  entire  range  for 
6e=6r  =  0°,  and  all  equilibria  are  stable  (Figure  3.16).  Note  how  roll 
rate  (p)  is  the  only  variable  reasonably  sensitive  to  6a  changes;  thus, 
there  is  also  decoupling.  But  the  rudder,  in  this  as  well  as  most  other 
aircraft  F  flight  regimes,  exhibits  far  different  characteristics.  This 
is  further  exemplified  by  the  time  history  shown  in  Figure  3.17.  This 
figure  shows  the  clear  growth  to  a  high-amplitude  limit  cycle  as  6r  steps 
in  10°  increments  from  0°  to  -40°,  with  6a  =  6e  =  0°.  Clearly,  the  effect 
of  rudder  is  nonlinear  in  this  region,  and  there  is  a  high  degree  of 
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coupling,  as  both  Figure  3.15  and  3.17  show. 

From  Point  A  in  Figure  3.15,  5r=-6.4°,  a  surface  is  generated  along 
the  5e  direction,  and  the  results  are  shown  in  Figure  3.18*  This  figure 
represents  a  parallel  slice  of  the  surface  only  6.4"  removed  along  the 
<5r  direction  from  the  projection  shown  in  Figure  3.10,  yet  it  is  evident 
how  dramatic  the  difference  is,  both  qualitatively  and  quantitatively. 

A  third  section  of  the  equilibrium  surface,  for  6r =  -25" ,  is  presented  in 
Figure  3.19,  and  represents  the  projection  alona  6e  centered  on  Point  B 
in  Figure  3.15.  Again,  there  is  a  considerable  difference  in  the  shape 
of  the  surface.  Note  further  that  inspecting  the  p-plot  alone  may  cause 
one  to  suspect  a  possible  bifurcation  point;  that  is,  the  point  where  the 
branches  intersect,  at  6e  -4°.  This  is  not  the  case,  however,  because 
there  is  no  such  overlap  for  the  (a,r)  curves. 

From  Point  C  in  Figure  3.18,  the  6e  =  0°  point,  an  orthogonal  pro¬ 
jection  is  generated  along  the  third  control  direction,  6a,  and  Figure  3.20 
shows  the  results.  For  the  (6e,6r)-values  given,  i.e.,  (0,-6, 4°),  the 
6a =  0°  point  is  a  precariously  stable  one.  Only  small  variations  are 
needed  to  cause  either  a  jump  (6a < 0)  or  development  of  an  unstable  or 
limit  cycle  conditions  (6a  >0). 

Returning  to  Figure  3.16,  for  which  r=0,  some  other  projections 

of  the  equilibrium  surface  are  compared  in  Figures  3.21  and  3.22.  In 

Figure  3.21,  with  the  same  controls  as  Figure  3.20,  i.e.,  6r=-6.4°, 

*In  this  figure,  and  in  others  where  it  appears,  an  asterisk  indicates 
that  two  pairs  of  complex  eigenvalues  of  [9f/9x]  have  positive  real 
parts--an  LL-condi tion. 
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6e  =  0°,  another  branch  of  the  equilibrium  surface  is  shown.  Aileron  has  now 
lost  its  linear  influence  here,  and  there  is  significant  coupling  as  well  as 
sharp  lack  of  symmetry  in  p.  Figure  3.22,  for  which  6r=-25°,  shows  very 
similar  curves  as  does  Figure  3.16,  for  which  6r  =  0;  however,  the  equili¬ 
brium  values  have  increased,  and  more  significantly,  roll  rate  is  almost 
totally  insensitive  to  6a  commands.  The  aircraft  is  in  an  autorotational 
state  in  roll.  Figure  3.23  is  a  case  similar  to  Figure  3.22,  i.e.,  6r=-25°, 
except  that  now  Se =  10.3°  (pitch  down).  It  can  be  seen  that  higher  equi¬ 
librium  values  ensue  for  this  6e,  with  somewhat  less  stability.  (Only  the 
r-fia  plot  is  shown  since  others  are  quite  similar). 

A  further  effect  of  the  rudder  can  be  seen  by  comparing  Figure  3.12 
with  Figure  3.24.  In  both  of  these  figures,  6e  =  -ll°  and  6a  is  the  inde¬ 
pendent  variable.  For  6r=  13°  in  Figure  3.24,  there  are  no  stable  regions 
on  the  branch  and  there  is  no  symmetry.  Both  of  these  results  are  expected. 
Elevator  influence  on  the  rudder  controllability  may  be  seen  by  comparing 
Figure  3.25  (6a  =15°,  6e  =  0°)  with  Figure  3.26  (6a  =15°,  6e=7.3°)  and 
Figure  3.27  (6a  =15°,  fie  =-11°).  These  surfaces  again  give  evidence  of 
the  richness  of  nonlinearities  and  hysteresis-type  behavior,  with  several 
different  kinds  of  equilibria  over  the  control  regime.  The  aileron  effect 
on  rudder  may  also  be  noted,  by  comparing  Figure  3.25,  for  which  6e = 0°  and 
6a  =15°,  with  Figure  3.15,  for  which  6e  =  0°  and  6a  =  0°.  There  is  some 
change,  although  not  to  the  degree  shown  by  the  elevator,  revealed  further  by 
a  comparison  of  Figure  3.26  (6e  =  7.3°)  with  Figure  3.25  (6e=0°)  especially 
in  a  and  p.  To  say  that  6r  has  nonlinear  influence  on  p  is  an  understate¬ 
ment,  after  inspecting  Figure  3.26.  Figure  3.27,  where  6e=-ll°  and 
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6a  =15°,  is  most  unusual  because  the  curves  close  upon  themselves  as 
rudder  is  varied.  The  geometry  is  more  evident  by  comparing  Figure  3.27 
with  Figure  3.28.  In  the  latter,  6r  =  0°,  6a  =  15°  and  6e  varies,  so  that 
it  is  "orthogonal"  to  Figure  3.27.  The  line  E-E'  represents  the  plane 
depicted  by  Figure  3.27;  conversely,  the  line  D-D'  in  Figure  3.27  is  in 
the  plane  depicted  by  Figure  3.28. 

The  equilibrium  surfaces  presented  in  this  section  were  shown  to 
provide  a  feel  for  the  great  variety  of  behavior  which  is  possible  in 
these  high-a  regimes.  One  can  begin  to  understand,  with  particular  ref¬ 
erence  to  the  time  history  results  shown  in  Fig.  3.29,  and  its  relevant 
equilibrium  surfaces  shown  in  Figs.  3.24  and  3.26,  that  in  certain  regions 
the  smallest  change  in  starting  conditions  can  result  in  widely  divergent 
results. 

3.3.3  Ming  Rock  Motions 

Wing  rock  has  been  mentioned  briefly  in  discussing  the  aircraft  F 
time  history  presented  in  Figure  3.11,  Section  3.3.1.  This  particular 
phenomenon  arises  as  the  result  of  developing  instability  of  airflow 
over  the  wings,  a  consequence  of  a  approaching  its  stall  value.  The 
main  feature  of  wing  rock  is  pronounced  roll  oscillations  whose  ampli¬ 
tude  increases  at  least  through  stall,  and  which  usually  couple  into 
yaw  and  pitch  oscillations.  The  coupling  effects  are  due  to  the  high-a 
nature  of  the  motion.  In  fact,  the  rolling  effects  so  prominent  in  wing 
rock  are  the  result  primarily  of  the  elevator  control  actions.  This  basic 
longitudinal-lateral  coupling  is  a  main  feature  of  high-a  dynamics. 

Figure  3.10  demonstrates  this  coupling  of  5e  deflections  into  the  lateral 
dynamics  (p,r),  and  should  be  referred  to  in  examining  the  results  shown 
in  Figure  3.11. 
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Figure  3.11a  shows  a  pre-stall  buffeting  which  has  only  traces  of 
lateral  oscillation.  Aileron  and  rudder  are  fixed  in  the  neutral  position 
for  Figure  3.11,  and  the  initial  state  is  in  trim.  As  can  be  seen,  velocity 
is  allowed  to  vary  in  3.11a  and  3.11b — the  aircraft  is  effectively  in  a  free-fall 
state--and  this  causes  damping  in  the  more  prominent  longitudinal  oscilla¬ 
tions.  In  Figure  3.11b,  the  damping  is  enhanced  by  changing  the  elevator 
control  sequence.  Here,  6a  is  set  and  held  at  -15°  after  20  seconds. 

However,  if  a  thrust  schedule  is  introduced  which  maintains  V  constant, 
the  6e  control  sequence  which  was  used  in  Figure  3.11a  produces  very  dif¬ 
ferent  results,  as  can  be  seen  in  Figure  3.11c.  In  this  figure,  the  roll 
rate  oscillations  become  very  severe,  and  they  induce  strong  pitch  oscillations 
as  well.  This  pronounced  limit  cycle  behavior  is  predicted  by  Figure  3.10, 
which  was  generated  assuming  constant  velocity. 

3.3.4  Post-Stall  Gyrations 

It  is  seen  from  inspection  of  the  global  equilibrium  surfaces  (e.g.. 

Figure  3.9),  which  show  both  spin  and  non-spin  regimes,  that  there  is  a 
basic  barrier  between  these  two  regimes.  This  is  due  physically,  in  part, 
to  the  great  changes  which  occur  in  the  state  of  the  vehicle  as  it  under¬ 
goes  transition  from  trim  to  spin  conditions.  The  velocity  vector  changes 
approximately  90°,  from  roughly  horizontal  to  vertical  (and  down) ;  angle- 
of-attack  similarly  undergoes  very  large  changes.  Time  history  runs  to 
be  shown  (Section  3.3.5  and  3.4)  indicate  that  if  the  initial  conditions  are 
close  to  either  those  of  stable,  developed  flat  spin  or  of  trim--and  the 
controls  are  set  to  values  representative  of  these  two  condi tions— then 
the  ensuing  motions  will  stay  in  these  regions  and  are  stable.  The  spin 
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motion  shows  a  wel 1 -developed  spiral,  with  constant  vertical  velocity, 
pitch,  roll  and  angle-of-attack  (for  an  aircraft  in  "flat"  spin,  the  equi¬ 
librium  pitch  angle  magnitude  is  typically  no  larger  than  about  15°;  thus, 
the  aircraft  is  spiraling  down  "on  its  belly"  with  substantial  yaw  rate).  But 
as  Fig.  3.9  shows,  if  the  controls  move  from  trim  (Point  D)  towards  spin  set¬ 
tings  (Point  A),  then  the  motions  indicative  of  stable,  developed  flat  spin  may 
not  result  since  limit  cycle  regimes  with  large  domains  of  attraction  exist 
around  Point  B.  Physically,  it  is  known  that  once  a  exceeds  its  stall  value, 
the  aircraft  becomes  subject  to  violent  oscillatory  motions  indicative  of 
the  loss-of-lift  condition  attendant  with  flow  separation  at  high-a. 

These  motions  are  called  post-stall  gyrations,  and  if  there  is  any  kind 
of  equilibrium  associated  with  them,  it  is  most  certainly  not  a  stable 
equilibrium  point.  Once  the  controls  move  to  values  where  high-*  non- 
linearities  predominate,  Hopf  Bifurcations  are  seen  to  occur  and  only  limit 
cycle  equilibria  exist.  And,  as  stated  earlier,  to  move  from  a  limit 
cycle  solution  to  a  stable  equilibrium  point  typically  requires  special 
sequencing  of  control  changes.  At  least  this  is  the  case  with  the  air¬ 
craft  F  model  as  defined  and  described  in  this  report.  It  is  easy  to 
sequence  the  controls  so  that  aircraft  F  enters  an  inadvertent  spin-like 
condition;  however,  this  spin  is  usually  not  the  smooth,  flat  spin  located 
at  Point  A  in  Figure  3.9.  The  entered  spin  is  predominantly  oscillatory, 
and  a  steep  or  intermediate  spin  which  corresponds  to  a  location  not  on  an 
equilibrium  surface  since  they  are  all  unstable  in  this  region,  but  one  evolv¬ 
ing  from  the  intermediate  region  of  the  post-stall  gyration  governed  by  limit 
cycle  (L)  branches.  This  region  is  basically  located  between  the  spin  and  non-spin 


89 


regions  and  is  featured  by  a  family  of  limit  cycles.  It  is  not  clear 

whether  these  limit  cycles  should  be  designated  as  "spin  motions"  or 

high-a  post-stall  gyrations,  as  they  cover  a  broad  range  of  a  values  (see  Fig.  3.9b). 

A  spin  entry  run  is  presented  in  Figure  3.29.  In  this  run,  the 
controls  were  initially  at  trim;  then  Se  was  set  to  -10°  at  20  sec,  6r  to 
-29°  at  40  seconds  and  6a  to  15°  at  t  =  45  seconds.  The  results  are  consistent 
with  what  has  been  discussed  above.  The  pitch-up  action,  pulling  back 
on  the  stick, causes  stall;  thereupon  lateral  control  inputs  trigger  the 
gyrational  limit  cycle  behavior. 

A  similar  run  is  shown  in  Figure  3.30  .  Here,  the  trajectory  begins 
with  trim  conditions,  but  6a  is  stepped  to  15°  and  6e  to  -11°,  and  held 
(rudder  remains  momentarily  at  0°).  The  initial  conditions  correspond  to 
the  S-segment  of  the  equilibrium  surface  shown  in  Figure  3.27.  As  6r  is 
then  varied  in  steps  to  14°,  this  induces  the  oscillatory  behavior  shown 
in  Figure  3.30.  The  mean  value  of  a  indicates  that  an  intermediate  spin 
has  been  achieved  (r  is  also  high,  with  q  very  small,  proportionally).  A 
spiral  of  sorts  has  most  likely  developed,  as  the  number  of  turns  (the 
plot  variable  TURNS)  is  increasing  at  a  steady  rate  (small  oscillations 
superimposed),  and  XN0RTH  and  YEAST  are  approaching  a  steady-state  mean 
value,  with  oscillation.  Finally,  a  rough  comparison  of  mean  values  of 
r  and  to  indicates  that  most  of  the  angular  motion  is  in  yaw. 

3.3.5  Spin  Equilibrium  Surfaces 

A  more  complete  study  of  spin  equilibrium  surfaces  will  be  given  in 
Sec.  3.4,  so  many  of  the  relevant  equilibrium  surfaces  will  be  shown  there. 


The  figures  presented  in  this  section  have  the  feature  that  they  were 
generated  using  the  full  eighth-order  spin  dynamic  system,  in  which  gravity 
effects  are  included  (pitch-roll  coupling)  and  velocity  is  allowed  to  vary. 

Figure  3.31  shows  the  r  and  a  vs.  6a  surfaces  for  a  left  pro-spin  control 
setting,  6e=0°,  6r  =  28.3°  (the  combination  6a  <  0°  and  6r>0°  produces  nega¬ 
tive  yaw  rate).  The  shape  of  these  curves  is  seen  to  be  quite  similar  to  the 
surface  generated  in  the  right  pro-spin  control  region  (see  Figure  3.34). 

Later,  when  the  non-spin  system  is  used  to  qenerate  the  spin  equilibrium 
surfaces,  it  will  be  seen  that  this  shape  persists,  although  the  numerical 
results  differ.  Figure  3.32  shows  that  the  absolute  variation  of  V  is  small 
(about  6%),  so  that  there  is  justification  in  assuming  V  constant.  A  surface 
projected  onto  the  r-axis  is  presented  in  Figure  3.33,  in  which  6r  is  the  in¬ 
dependent  control  and  6a  =  -15°,  6e  =  0°  (left  pro-spin  controls).  Noting 
again  the  common  shape  vis-a-vis  Figure  3.31,  the  6a  plot,  and  the  fact 
that  the  right  pro-spin  control  region  possesses  the  same  type  of  mani¬ 
fold,  the  right  and  left  pro-spin  manifolds  are  presented  over  the  6a-6r 
plane  as  shown  in  Figure  3.34.  Manifold  A  represents  the  right  pro-spin 
manifold  (surface)  and  Manifold  B  represents  the  left  pro-spin  manifold.  Non-spin 
equilibrium  surfaces  are  not  shown  in  this  figure;  they  would  be  centered 
about  the  origin  and  would  not  be  in  contact  with  either  spin  manifold. 

While  aircraft  F  possesses  symmetry  to  the  extent  displayed  by  the 
presence  of  two  spin  manifolds  of  similar  shape,  it  will  be  noted  that  the 
left  oro-spin  branches  presented  in  this  section  do  not  possess  stable  segments, 
whereas  the  right  pro-spin  branches,  shown  in  Sec.  3.4,  do.  This  is  a 
reflection  of  the  asymmetry  of  the  aerodynamic  data  which  is  used  by  the  model. 


91 


3.3.6  Spin  Entry  Time  Histories 

Given  that  a  stable,  developed  spin  equilibrium  manifold  exists  (for 
right  pro-spin  controls),  the  problem  of  reaching  this  manifold  from 
non-spin  flight  conditions  remains  to  be  considered.  As  discussed  earlier, 
there  exists  a  large  intermediate  region  between  the  clearly-defined  spin 
region  and  the  non-spin  region.  It  corresponds  to  the  flight  regime  which 
is  often  categorized  "post-stall  gyrations."  The  aircraft  motion  in 
this  regime  sometimes  appears  chaotic,  with  large  oscillations  of  often- 
irregular  frequency  but  is  mostly  of  a  limit  cycle  type.  We  will  charac¬ 
terize  this  region  as  a  limit  cycle  region  based  on  such  results  as  pre¬ 
sented  in  Figure  3.9.  It  effectively  acts  as  a  barrier  to  a  sudden  jump 
from  trim  to  spin  conditions  or  vice  versa.  The  time  histories  presented 
in  this  section  show  this  behavior,  and  they  further  indicate  that  a  large 
segment  of  this  intermediate,  limit  cycle  region  is  characterized  by 
motions  normally  designated  as  oscillatory  spin.  Furthermore,  it  is 
usually  steep  (a*  55°)  or  intermediate  (a =70°)  in  nature,  based  on  the 
mean  value  of  angle-of-attack. 

Figure  3.35  presents  a  case  in  which  the  controls  are  moved  from  neutral 
to  §5pjN  at  t=2  seconds.  (The  initial  flight  condition  is  trim.  £SPjN 
is  defined  at  the  beginning  of  Sec.  3.3.)  If  the  spin  manifolds  were  as 
simple,  relatively,  as  those  of  the  roll  departure  region,  say,  then  one  would 
expect  to  see  entry  into  developed  spin.  This  does  not  happen,  however.  There 
is  clearly  a  post-stall  condition,  but  yaw  rate  does  not  reach  the  required  level 
(about  lOOVsec),  and  generally  the  energy  interchange  maintains  a  somewhat  chaotic 
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post-stall  gyration.  In  this  run,  velocity  was  allowed  to  vary.  If,  however, 

V  is  fixed  to  the  value  corresponding  to  the  stable  segment  of  the  equilibrium 
spin  manifold  for  the  control  set  6spiN  (V  =  443  fps),  then  the  rerun  of  the 
above  trajectory  produces  a  more  uniform  transient  oscillatory  behavior  which 
decays  somewhat  towards  a  mildly  oscillatory  spin  condition  (Figure  3.36  ). 

This  example  points  out  that  the  role  of  thrust  in  spin  entry  studies  must 
receive  more  attention.  By  forcing  V  =  0  we  are  effectively  maintaining 
thrust  at  a  level  which  exactly  opposes  aerodynamic  drag.  This  becomes 
physically  unrealistic,  however,  in  simulations  where  a  and  B  undergo  dramatic 
variations,  as  happens  here;  this  means  that  the  thrust  direction  is  fluctuating 
wildly,  as  well  as  thrust  magnitude. 

The  maneuver  discussed  with  regards  to  Fig.  3.12  and  3.13  may  also 
be  considered  as  a  spin  entry  sequence.  Fig.  3.37  shows  a  variation  of 
this  maneuver  in  which  the  rudder  takes  the  place  of  the  aileron  as  the 
lateral  control  during  the  maneuver.  As  previously  mentioned,  the  rudder 
plays  a  more  critical  role  in  converting  post-stall  motion  into  spin 
entry  and  subsequent  spin  motion,  and  Fig.  3.37  shows  that  an  oscillatory, 
left  pro-spin  is  induced  by  maintaining  <5r  at  10°  while  the  elevator  is 
stepped  from  0°  to  -14°  in  2  second  intervals,  2°  at  a  time  from  t=0.  It 
is  not  in  this  case  necessary  to  step  5e  in  sequences  to  -14°  to  show  this 
effect,  but  we  have  done  so  here  and  in  other  instances  to  observe  the 
effect  of  intermediate  control  values.  However,  more  care  is  needed  in 
general  with  large  control  changes,  as  unwanted  jump  phenomena  may  occur, 
placing  the  motion  in  a  different  region  than  desired. 

The  following  runs  present  variations  on  the  run  depicted  in  Fig.  3.36. 

In  these  runs,  the  control  sequencing  for  spin  entry  was  made  more  realistic: 
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5e  =  -21°  from  t  =  0  sec;  <5r  changed  from  0°  to  -25°  at  t  =  2  sec;  and  6a 

changed  from  0°  to  15°  at  t =  6  sec.  In  addition  V  was  fixed  at  443  fps  for 

all  runs  and  the  trajectory  begins  from  trim  conditions,  as  does  the  case 
shown  in  Fig.  3.36.  The  run  shown  in  Fig.  3.38  is  exactly  similar 
to  the  Figure  3.36  case,  except  for  the  difference  in  control  sequencing 
mentioned  above.  Comparison  shows  that  the  final,  oscillatory  state  is 
quite  similar. 

Following  the  example  of  Bihrle  (1976),  the  next  runs  show  the  effect 
of  changing  the  initial  roll  angle  from  a  trim  value  to  a  number  such  as 
60°.  Figure  3.39  shows  that,  again,  the  only  significant  difference  is  in 
the  transient  region,  which  lasts  about  20  seconds.  If  initial  pitch  angle 

is  changed  to  -50°  in  addition.  Figure  3.40,  the  ensuing  motion  is  sub¬ 

stantially  different.  Yaw  rate  does  not  achieve  the  same  value,  and  neither 
does  angle-of-attack.  These  results  indicate  that  more  study  of  the 
effects  of  changing  the  initial  state  would  be  desirable.  As  an  eventual 
goal,  methods  for  determining  the  domains  of  attraction  of  all  stable 
equilibria  (points  and  limit  cycles)  should  be  developed. 

3.4  Developed  Spin  Motion 

It  has  been  mentioned  earlier  that  the  task  of  effecting  transition 
from  non-spin  to  stable,  flat  spin  equilibrium  requires  passing  through 
a  highly  chaotic  and  oscillatory  intermediate  region  of  post-stall  gy¬ 
rations.  It  is  consequently  much  easier  to  study  spin  behavior  by 
making  time  history  runs  whose  initial  conditions  correspond  to 
^-SPIN’-SPIN^'  tlie  trajectory  begins  at 


I 


94 


-SPIN  ( P  »Q » i"  tO.  ,  B , V ,  0  ,<}> ) 

=  (30,-4,100,73.5,-3,443,-16.6,-2.29) 
and 

§  =  6srin  *  (<5a,Se,6r)  =  (15,-21,-25) 

(all  angular  terms  in  degrees,  V  in  fps),  the  ensuing  spin  motion  is  very 
smooth,  as  indicated  by  Figure  3.41.  This  figure  shows  the  horizontal 
trace  of  the  vehicle  center  of  mass,  and  the  altitude  variation  for  the 
initial  conditions  described  above.  These  conditions  are  in  the  middle 
of  the  stable  spin  equilibrium  branch.  Point  A  of  Figure  3.42.  At  t  =  15 
seconds,  6r  is  changed  to  -29°,  so  that  a  jump  occurs  to  the  upper  limit 
cycle  branch,  Point  B.  The  only  apparent  result  of  this  6r  change  in  Fig¬ 
ure  3.41  is  a  slight  tightening  of  the  spiral,  although  yaw  rate  increases 
more  dramatically.  Figure  3.43,  which  time-shifts  the  <5r  change  to  -29°  at 
t=0,  does  confirm  the  growth  to  a  limit  cycle  condition. 

From  Figure  3.41  it  can  be  observed  that  the  equilibrium  flat  spin  for 
aircraft  F  generates  a  very  tight  spiral  which  slowly  drifts  to  the  right, 
due  to  asymmetries  in  the  aerodynamic  data.  Also,  it  can  be  seen  that  the 
rate  of  descent  is  constant.  Additionally,  although  it  is  not  shown  here, 
the  state  variables  remain  constant  for  <5r  =  -25°,  and  exhibit  well-damped 
transient  behavior  to  new  steady  state  values  when  6r  is  changed  to  -29°. 
Virtually  all  of  the  vehicle's  velocity  in  this  condition  is  vertical  and 
the  aircraft's  orientation  with  respect  to  the  horizon  is  almost  "flat" 

(9=  -16°,  <p  =  -3°),  with  a  very  high  yaw  rate  (about  100  degrees  per  second). 
The  following  sections  discuss  these  and  similar  results  in  more  detail. 
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3.4.1  Equilibrium  Surfaces  in  the  Developed  Spin  Region 

The  typical  spin  equilibrium  curve  is  seen  in  Figure  3.42  .  In  this 
figure,  Sr  is  varied  while  Sa  =  15°  and  6e  =  -21°  and  V  is  fixed  at  443  fps. 

A  very  similar  curve  is  shown  in  Figure  3.44  .  Here,  6e  =  0°  but  that  is 
the  only  difference  from  Figure  3.42  .  Obviously,  the  elevator  has  little 
effect  on  this  particular  type  of  spin  motion.  Both  of  these  surfaces  were 
generated  using  the  non-spin  equilibrium  set,  neglecting  gravity. 

Choosing  6r=  -25°,  the  projection  along  6a  is  shown  for  a  in  Figure  3.45. 
Again,  V  =  443,  6e  =  0°  and  the  non-spin  set  (which  neglects  gravity)  is 
used.  If  6e  is  changed  to  -21°,  the  a  vs.  6a  branch  shown  in  Figure  3.46 
results. 

Finally,  a  composite  of  all  of  the  relevant  equilibrium  points,  for  all 
relevant  control  states,  is  projected  onto  the  (r-ct)  plane  (Fig.  3.47).  These 
variables  are  the  most  significant  ones  in  terms  of  analyzing  spin  motion. 

Very  noticeable  is  the  "gap"  between  the  non-spin  and  spin  regions. 

3.4.2  Importance  of  Assumptions  Concerning  Spin  Equilibria;  Comparisons 

In  this  section  the  significance  and  validity  of  some  of  the  assumptions 
relevant  to  generating  the  spin  region  equilibrium  surfaces  is  discussed. 

This  is  done  mainly  by  means  of  comparison  of  various  effects. 

Une  result  mentioned  in  the  last  section  is  that,  once  the  aircraft 
is  in  the  stable  flat  spin  condition,  elevator  controllability  becomes 
negligible.  This  is  readily  seen  by  comparing  Figure  3.42,  where 
6e  =  -21°,  with  Figure  3.44  where  6e  =  0°.  Figure  3.48  shows  the 


small  effect  of  6e  in  this  spin  equilibrium  region.  Only  roll  rate  is 
moderately  affected  as  6e  takes  values  of  0°,  -11°  and  -21°;  and  angle- 
of-attack,  surprisingly,  is  effectively  unchanged.  A  possible  explanation 
is  the  large  value  of  local  sideslip  at  the  elevator  locations,  generated 
by  the  large  steady  yaw  rate.  Notice  that  V = 600  in  this  figure. 

Of  somewhat  more  importance,  however,  is  the  validity  of  the  assumption 
of  fixing  velocity  and  neglecting  gravity  effects  in  the  spin  region.  If 
this  assumption  can  be  accepted  as  valid  for  initial  phases  of  analysis  of 
spin  motion,  then  spin  region  equilibrium  and  bifurcation  surfaces  can  be 
generated  by  the  simpler  5  DOF  non-spin  system  of  equations.  Figure  3.49 
shows  that,  at  least  for  aircraft  F  in  the  flat  spin  region, 
the  simplifying  assumption  V  -  const,  g  =  0  may  be  used  for  initial  spin 
analyses.  In  fact,  changing  the  velocity  is  seen  to  produce  greater  dif¬ 
ferences  than  neglecting  gravity.  In  this  figure,  6a  =  15°  and  6e  =  -21°; 
also,  the  g  f  0  branch,  because  it  was  run  using  the  full  spin  system, 
does  not  have  associated  with  it  a  constant  velocity.  However,  as 
Figure  3.32  shows,  the  V  range  is  only  about  20  fps.  Another  surprising 
result  is  the  total  insensitivity  of  angle-of-attack  to  these  changes. 

A  conclusion  to  be  drawn  from  this  comparison  is  that,  for  flat 
developed  spin  using  the  aircraft  F  model,  since  pitch  angle  is  small 
(about  15°)  and  a  =  90°,  the  term  cos  a  sin  6,  which  couples  gravity  into 
the  a-equation,  is  quite  small.  Furthermore,  gravity  does  not  directly 
couple  into  the  p  and  r  equations,  as  does  V  (through  dynamic  pressure  and 
by  influencing  the  aerodynamic  coefficients). 


97 


Another  observation  of  practical  significance  is  that  only  one  consis 
tent  set  of  assumptions  is  needed  in  order  to  generate  a  truly  global 
equilibrium  surface,  examples  of  which  are  shown  in  Figs.  3.9,  3.50  and 
3.51.  One  needs  only  to  choose  a  "reasonable"  value  for  V  in  order  to 
generate  these  figures  with  the  easier  non-spin  equations.  Notice  that 
it  is  not  possible--and  we  did  make  simulation  runs  to  verify  this--to 
use  the  spin  equations,  with  nonzero  gravity,  to  produce  equilibrium  mani¬ 
folds  in  non-spin  regions;  roll  and  pitch  angles,  which  must  be  included 
in  this  system,  have  no  equilibrium  values  in  most  non-spin  regions. 
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3. 5  Spin  Recovery  and  Prevention 

In  terms  of  the  concepts  employed  by  BACTM  to  analyze  aircraft 
behavior,  it  is  possible  to  state  the  goal  of  spin  recovery  as  follows: 
spin  recovery  is  achieved  by  control  sequences  which  move  the  equilibrium 
point  from  a  stable  spin  region  to  either  an  unstable  spin  equilibrium 
point  or  a  point  on  a  non-spin  equilibrium  branch,  stable  or  unstable,  from 
which  other  control  actions  can  produce  trim  conditions.  A  jump  from  one 
stable  equilibrium  to  another  in  the  spin  regime  is  undesirable. 

With  regard  to  prevention  of  spin  situations,  it  will  be  seen  that 
the  rudder  is  the  most  sensitive  aerosurface  control  for  aircraft  F,  in 
terms  of  spin  entry.  The  aileron  also  has  considerable  influence,  but 
with  this  particular  model,  rudder  influence  predominates.  In  view  of 
this,  high-speed/large-attitude-change  maneuvers  which  require  large 
elevator  and  aileron  deflections  (e.g.,  a  rolling  pull-up  maneuver)  be¬ 
come  very  prone  to  spin-entry  situations  unless  use  of  the  rudder  is 
carefully  controlled. 

The  next  section  deals  with  some  of  the  aspects  of  spin  recovery, 
based  on  BACTM  analysis  using  aircraft  F,  and  the  following  section  will 
cover  aspects  of  spin  prevention. 

3.5.1  Spin  Recovery  with  Aircraft  F 

For  aircraft  F,  a  right  pro-spin  control  setting  of  (6a,6e,6r)=(15°,-21°,-25°) 
designated  §spj ^ *  represents  the  "spin  control  setting."  This  setting, 
along  with  the  proper  values  of  the  state  variables,  x^pjN»  results  in  a 
flat,  equilibrium  spin.  See  Point  A  in  Fig.  3.42.  The  magnitude  and  sign 
sense  of  the  controls  §SPIN  is  very  representative  of  similar  spin  settings 
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of  other  aircraft;  the  elevator  (6e)  is  large  and  negative,  to  provide  the 
high  a  needed  for  stall  and  subsequent  spin  (although  once  in  stable  spin, 
the  aircraft  F  equilibrium  state  is  quite  insensitive  to  elevator  control 
actions--see  Figure  3.48);  the  aileron  (6a)  is  at  its  extreme  setting, 
of  opposite  sign  to  the  rudder  (6r);  and  the  large  negative  rudder 
generates  the  high  positive  yaw  rate  which  signifies  the  development 
of  spin  behavior  following  post-stall  gyrations.  The  aileron  is  of 
opposite  sign  in  a  spin  setting  because  of  the  effect  of  adverse  yaw 
due  to  the  aileron;  that  is,  for  a  positive  yaw  rate  ( 6 r  <  0 ) ,  a  negative 
roll  rate  (6a  >0)  induces  a  positive  yaw  moment,  thereby  enhancing  the  yaw 
rate.  In  a  coordinated  turn,  both  6a  and  Sr  have  the  same  sign,  and  yaw  rate 
(r) is  predominantly  sensitive  only  to  5r.  Proper  sequencing  of  controls  for 
spin  entry  is  important,  because  hysteresis  effects  are  especially  pro¬ 
minent  in  these  high-t  regions.  For  example,  the  elevator  is  a  much  more 
effective  control  for  spin  entry  when  it  is  applied  while  sideslip  (6)  is 
still  small  in  magnitude.  This  is  consistent  with  the  usual  circumstance 
of  spin  following  a  stall;  and  the  elevator  typically  induces  stall  be¬ 
cause  of  its  direct  influence  on  angl e-of-attack.  See  Fig.  3.10.  Further¬ 
more,  elevator  control  effectiveness  is  practically  nonexistent  for  high 
values  of  6. 

When  a  significant  yaw  rate  is  added  to  post-stall  motions,  and  when 
angl  e-of-attack  (a)  approaches  values  for  which  autorotation  in  yaw  is 
possible  (i.e.,  Cn~0),  then  entry  into  spin  will  most  likely  result.  For 
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aircraft  F,  Cn  =  0  around  a  =  65°.  Spin  may  be  considered  a  form  of  auto- 
rotational  yawing.  This  is  a  condition  marked  by  a  general  ineffective¬ 
ness  of  the  lateral  controls.  However,  by  using  the  equilibrium  surfaces, 
recovery  control  sequences  can  be  developed. 

The  standard  method  for  spin  recovery  is  to  rapidly  proceed  to  an  anti¬ 
spin  setting--i.e. ,  for  our  example,  this  would  involve  zeroing  the  elevator 
and  fixing  6r  and  6a  at  their  opposite  extremes.  This  and  similar  techniques 
based  on  aerosurface  control  actions  alone  are  not  always  effective,  and  air¬ 
craft  are  often  equipped  with  special  thrusters  and  drag  parachutes  for  spin 
recovery  purposes.  However,  it  is  possible  to  effect  recovery  from  spin  with 
the  aircraft  F  model,  and  Fig.  3.9  indicates  how  this  may  be  done.*  This  fig¬ 
ure  shows  the  equilibrium  surfaces  for  6a  =  15°,  6e  =  0°,  V=  600  and  g  =  0. 

From  the  spin  state,  the  elevator  is  returned  to  the  neutral  position; 

this  corresponds  to  point  A  in  Fig.  3.9  .  Then,  the  rudder  is  increased 

from  -25°  to  at  least  15°.  This  will  induce  two  jumps,  the  first  one  to 

a  limit  cycle  around  point  B,  and  the  second  to  a  limit  cycle  around 

point  C  on  the  lower  branch.  Then  6r  may  be  decreased  to  0°  (point  D), 

as  this  last  equilibrium  branch  passes  through  small  values  for  a  and  r. 

These  control  actions  must  be  taken  over  a  long  enough  period  to  allow  the 

transient  motions  to  die  out.  Point  D  in  Fig.  3.9  corresponds  to  point  D  in 

Fig.  3.16,  which  is  an  equilibrium  surface  showing  the  final  recovery  sequence: 

roll  rate  p  is  reduced  by  returning  the  aileron  to  its  neutral  setting.  Note  that  it 

involves  a  control  effort  (change  in  rudder  from  one  extreme  setting  to  the 

other)  which  is  used  in  practice  for  spin  recovery.  Again,  the  equilibrium 

surfaces  generated  by  BACTM  tend  to  confirm  previous  results  and  past 

♦There  are  other  possibilities  which  we  shall  try  to  investigate  at  a  later 
time,  with  more  comprehensive  models,  such  as  the  F-4. 
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experience,  while  at  the  same  time  pointing  out  alternative  possibilities 
(which  we  expect  to  explore  in  more  detail  at  a  later  time).  A  further 
characterization  of  the  limit  cycle  behavior  is  required  for  obtaininq  op¬ 
timal  recovery  procedures.  This  is  because  more  information  needed  about 
the  limit  cycle  domains  of  attraction  in  the  sensitive  intermediate-a 
region  (25°  s a  £  65°).  Eventually,  a  complete  calculation  of  the  bifurca¬ 
tion  surfaces  (both  Hopf  and  elementary)  should  be  done  for  spin  recovery 
control  design. 

Another  possibility  for  a  spin  recovery  strategy  "leading  with  5a"  is 
presented  in  Fig.  3.50.  Here,  6a  is  reduced  from  15°  so  that  a  jump  occurs, 
to  a  limit  cycle  condition  around  point  B  from  point  A  (the  latter  point  is 
the  same  point  A  in  Fig.  3.9).  This  jump  will  actually  increase  p,  the  roll 
rate,  but  this  is  a  desirable  method  of  rolling  the  aircraft  into  the 
airflow,  which  reduces  a,  as  can  be  seen.  The  final  step,  then,  is  to 
reduce  6r  in  magnitude  to  its  neutral  setting,  and  Fig.  3.15  shows  that 
this  returns  (r,a,p)  to  trim  values.  Similarly,  time  history  runs  confirm 
this.  As  with  Fig.  3.9,  Fig.  3.50  represents  equilibrium  conditions 
for  6e  =  0°,  V  =  600  fps  and  g  =  0. 

It  should  be  recognized  that  excessive  application  of  the  "anti-spin" 
control  actions,  in  attempting  to  effect  spin  recovery,  can  lead  to  a 
"reverse  spin"  situation  if  the  controls  are  not  moved  towards  neutral 
quickly  enough.  A  common  method  of  spin  recovery  is  to  oscillate  the 
controls,  particularly  6a  and  6r,  back  and  forth  between  their  limits.  The 
frequency  of  the  oscillations  is  usually  determined  by  visual  cues  avail¬ 
able  to  the  pilot,  in  conjunction  with  the  autopilot.  Fig.  3.51  shows  a 
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spin  reversal  situation.  At  42  seconds,  the  controls  were  changed  from 
right  pro-spin  (6a,Se,6r) =  (15°, -21°, -30°),  to  left  pro-spin,  (-15°, -21°, 
It  can  be  seen  that,  within  5  or  6  seconds,  the  yaw  rate  has  changed 
sign  but  not  magnitude,  and  the  basic  spin  condition  remains  otherwise 
unaffected. 
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3.6  Explanation  of  Spin  Behavior  of  Aircraft  F 

As  with  other  aircraft  analyzed  in  the  literature  for  spin  behavior, 
the  main  feature  of  aircraft  F  in  developed  spin  is  an  extremely  high 
angle-of-attack  and  persistent,  steady  yaw  rate.  The  presence  of  these 
conditions  simultaneously,  without  major  fluctuation  between  high  and  low 
values  of  r  and  a,  indicates  the  spin  condition.  Once  the  aircraft  has 
been  maneuvered  into  a  stall  condition,  both  the  equilibrium  surfaces  and 
the  time  history  simulations  indicate  that  wild,  oscillatory  post-stall 
oscillations  and  gyrations  result.  The  aircraft  has  entered  a  flight 
regime  lacking  in  stable  equilibrium  points;  if  the  lateral  controls  have 
been  set  to  "pro-spin"  positions  just  prior  to  or  at  the  onset  of  stall, 
then  the  yawing  motion  will  predominate  the  post-stall  dynamics.  If  inertia 
coupling  and  aerodynamic  forces  and  moments  are  then  phased  together  so 
that  the  values  of  angle-of-attack  and  sideslip  generate  negligible  yaw- 
moment  coefficient,  Cp,  then  the  yawing  motion  will  become  autorotational . 
For  aircraft  F,  this  will  occur  when  a  is  about  70°  and  0  is  within  ±10°. 

In  a  "pro-spin"  control  setting,  the  aileron  is  moved  in  the  opposite 
sense  to  the  rudder,  and  both  controls  are  typically  at  or  close  to  their 
extreme  values.  Thus  if  aileron  is  positive,  for  negative  roll  rate,  but 
the  rudder  is  negative,  the  positive  yaw  rate  generated  by  the  rudder 
will  be  enhanced  by  the  adverse  aileron  yaw. 

The  transition  dynamics  from  trim  condition  to  spin  equilibrium  for 
aircraft  F  involves  limit  cycle  oscillations.  A  thorough  analysis  of  the 
generation  and  interplay  of  the  aerodynamic,  inertial  and  gravitational 


forces  and  moments  is  not  merited  for  this  model,  since  aerodynamic  data  were 
developed  in  order  to  simulate  full  spin  motion.  However,  it  is  felt  that 
post  stall  gyrations  produced  by  the  model  are  similar,  in  general, 
to  what  is  actually  experienced  in  flight  tests  of  military  aircraft. 

The  model  limitations  and  open-loop  nature  of  the  simulation  make  complete 
transition  from  trim  conditions,  controls  neutral,  to  stable,  developed  spin 
a  very  difficult  task.  The  equilibrium  surfaces  are  useful  starting  points, 
but  in  regions  where  stable  branches  are  nonexistent,  they  cannot  predict 
easily  the  nature  of  the  motion  to  be  encountered  in  that  region.  Thus,  an 
exhaustive  series  of  runs  would  be  required  in  order  to  proceed  completely 
from  trim  to  spin.  However,  we  have  been  able  to  show  that  by  making  two  runs, 
one  with  initial  conditions  at  UsPIN’-SPIN^  anc*  the  other  at  trim,  that 
respective  control  sequences  may  be  chosen,  using  the  information  provided 
by  the  equilibrium  surfaces,  so  as  to  bring  them  to  a  common  final  state. 

Fig.  3.52  shows  a  trajectory  which  began  in  a  spin  condition.  At  0.5  seconds, 
the  rudder  was  changed  from  -25°  to  25°,  and  6a  changed  from  15°  initially  to 
7.5°  (t  =  20  seconds)  and  0°  at  30  seconds.  The  resulting  oscillatory  spin 
matches  very  closely  the  final  condition  of  the  trajectory  shown  in  Fig.  3.53, 
which  began  in  trim  and  had  the  following  control  sequence: 
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This  is  an  important  result  in  that  it  does  show  that  the  model  produces 
motions  that  proceed  from  trim  to  spin  and  vice  versa.  As  the  equilibrium 
surfaces  show,  there  exist  several  attractors  (both  points  and  orbits)  and 
their  associated  domains  of  attraction  which  make  it  very  difficult  to 
effect  excursions  in  the  state-control  space.  We  have  effectively  com¬ 
bined  two  runs  which  take  advantage  of  more  favorable  equilibrium  struc¬ 
tures  in  their  respective  starting  regions  and  brought  them  to  a  common 
point.  From  this  point,  completely  new  control  sequences  must  be  used  in 
order  to  return  to  either  starting  point.  The  composite  equilibrium  sur¬ 
faces  shown  in  Fig.  3.9  and  3.50,  for  6e  =  0°,  show  that  the  jump  form  the 
spin  condition  to  highly  oscillatory  regions  is  much  easier  than  going 
in  the  reverse  direction.  Here,  it  can  be  seen  that  the  jump  is  to  a 
limit  cycle  surface  in  the  (6a, 6r)  plane.  Once  the  jump  occurs,  elevator 
controllability  becomes  more  prominent,  as  can  be  seen  by  comparing  Fig.  3.9 
with  Fig.  3.54,  which  has  6e=-ll°.  The  elevator  change  does  not  appre¬ 
ciably  affect  the  shape  of  the  spin  equilibrium  branch,  but  greatly  changes 
the  non-spin  branch. 

The  ease  with  which  one  can  move  from  one  point  to  the  other  along 
the  equilibrium  surfaces  is  influenced  greatly,  as  mentioned  above,  by  the 
structural  stability  properties  of  these  surfaces  in  control  space.  An¬ 
other  related  factor  is  the  location  of  the  attractors  in  this  space  and 
the  various  "domains  of  attraction."  Based  on  our  results,  it  seems  that 
the  domain  of  attraction  for  the  stable  segment  of  the  spin  equilibrium 
branch  is  much  smaller  than  neighboring  domains  of  attraction  of  limit 
cycles.  This  necessitates  precision  control  sequencing,  sensitively 
attuned  to  the  current  state,  in  order  to  achieve  stable  spin. 
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FIGURE  3.9(a) 

Equilibrium  Surface:  r,cx,p  vs.  6r 
5a  =  15°,  <5e  =  0° ,  V=  600  fps,  g=0 
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FIGURE  3.9(b)  (cont.) 
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FIGURE  3.10(a) 

Equilibrium  Surface:  r,o,p  vs.  fie 
fia  =  0°,  fir  =  0° ,  V  =  600  fps,  g  =  0 

The  *  refers  to  an  LL  instability 
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FIGURE  3.10(b)  (cont.) 
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FIGURE  3.10(c)  (concluded) 
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6a,  6r,  6e  same  as  3.11(a);  V=600  fps 


FIGURE  3.12(a) 

Equilibrium  Surface:  r,a,p  vs.  6a 
Se  =  -11°,  6r  =  0° ,  V  =  600  fps,  g  =  0 
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FIGURE  3.12(b)  (cont.) 
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FIGURE  3.12(c)  (concluded) 
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FIGURE  3.14(a) 

Equilibrium  Surface:  r.a.p  vs.  Sr 
sa  =  0°,  6e  =  -20° ,  V  =  600  fps,  g  =  0 
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FIGURE  3.15(a) 

Equilibrium  Surface:  r,a,p  vs.  6r 
6a  =  0°,  6e  =  0° ,  V  =  600  fps,  g  =  0 
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FIGURE  3.15(b)  (cont.) 
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Figure  3.18(a) 

Equilibrium  Surface:  r,a,p  vs.  6e 
V=  600  fps,  g=  0. 

The  symbol  *  refers  to  an  LL-type  instability  (two  unstable  complex  pairs) 
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Figure  3.18(b)  (cont.) 


Dflr 

DR= 


0-0 

-6.4 


140 


S  D  UL'UUL)  J  LUUU  Q  LLL'JU 

U  E  UUUULIU  K  LUUUU  R  LLLIUUU 

L  F  UUUUUUU  M  LUUUUU  X  LLlU 

h  UU  G  UUUUUUUU  N  LUUUUUU  v  llluu 

S  IJU'J  H  LU  0  LLU  Z  LLi.L 

g  C  UUUU  [  LUU  P  LLUU 

C. 

CD  I 


I 


?4 — 

-30.30 


-2?Too  ^M.no  -p.on 

DE 


- 

?.oo 


4  -- — - 1 

ip.od  ip.-x 


Figure  3.18(c)  (concluded) 
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FIGURE  3.19(c)  (concluded) 
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Figure  3.22(a) 

Equil ibrium  Surface:  r,a,p  vs.  6a 
V  =  600  fps,  g  =  0 
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Figure  3.22(b)  (cont.) 
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Figure  3.22(c)  (concluded) 
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FIGURE  3.23 

Equilibrium  Surface:  r  vs.  6a 
5e  =  10.3°,  6r  =  -25° ,  V  =  600  fps,  g  =  0 
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FIGURE  3.25(a) 

Equilibrium  Surface:  r,a,p  vs.  <5r 
6a  =  15°,  6e  =  0°,  V  =  600  fps,  g  =  0 
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FIGURE  3.25(c)  (concluded) 
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FIGURE  3.26(a) 

Equilibrium  Surface:  r,a,p  vs.  6r 
6a  =  15°,  6e  =  7.3°,  V  =  600  fps,  g  =  0 
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FIGURE  3.26(c)  (concluded) 
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FIGURE  3.27(a) 

Equilibrium  Surface:  r,a,p  vs.  6r 
6a  =  15°,  6e  =  -11° ,  V  =  600  fps,  g  =  0 
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FIGURE  3.27(b)  (cont.) 
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FIGURE  3.28(a) 

Equilibrium  Surface:  r,a,p  vs.  Se 
<5a  =  15°  ,  6r  =  0° ,  V  =  600  fps ,  g  =  0 
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Figure  3.?9 

Time  History:  Spin  Entry  Study 
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Figure  3.31(b)  (concluded) 

Equilibrium  Surface  -  Spin  Regime:  Angle  of  Attack  vs.  Aileron  Deflection 
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Figure  3.32 

Equilibrium  Surface  -  Spin  Regime:  Velocity  vs.  Aileron  Deflection 
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Figure  3.33 

Equilibrium  Surface  -  Spin  Regime;  Yaw  Rate  vs.  Rudder  Defection 
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Figure  3.34 

Equilibrium  Manifolds  -  Spin  Regime,  Aircraft  F 
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Figure  3.38(a) 

Time  History:  Spin  Entry;  6e=-21°,  V  =  0,  V=  443  fps 
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Figure  3.39(c)  (cont.) 
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Figure  3.44(c)  (concluded) 
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FIGURE  3.45 

Equilibrium  Surface:  a  vs.  6a 
6e  =  0°,  6r= -25°,  V  =  443  fps,  g  =  0 
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Figure  3.47 

Equilibrium  Surfaces  -  Composite  Angle  of  Attack  vs.  Yaw  Rate  at 
Various  Control  Settings 
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Equilibrium  Surface:  Comparison  of  Elevator 
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FIGURE  3.48(c)  (concluded) 
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FIGURE  3.49(a) 

Equilibrium  Surfaces:  Comparison  of  V,  g  Effects 
r,a,p  vs.  6r;  6a  =  15°,  6e  =  -21° 
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FIGURE  3.49(b)  (cont.) 
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FIGURE  3.49(c)  (concluded) 
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Figure  3.52(b)  (cont. ) 
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FIGURE  3.53(b)  (cont.) 
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CHAPTER  IV 


Other  Topics 

During  this  reporting  period,  further  work  has  been  performed  to 
understand  the  nature  of  the  nonlinear  high-a  dynamic  behavior  of  air¬ 
craft  H  (Mehra  et  al.  (1977)).  In  this  regard,  there  are  two  topics 
which  merit  discussion  at  this  time.  The  first  topic,  exploring  the  power 
spectra  of  some  motions,  is  part  of  an  overall  study  which  will  be  ex¬ 
panded  on  in  the  future  with  the  F-4  aircraft  model.  The  second  topic 
deals  with  using  our  knowledge  of  the  global  (nonlinear)  characteristics 
of  the  ari craft  model  of  interest  to  synthesize  a  command/stability  aug¬ 
mentation  system. 

4.1  Power  Spectra  of  Time  Histories  for  A/C  H 

As  a  means  of  qaining  further  insight  into  the  nature  of  the  Hopf 
Bifurcation  and  the  limit  cycle  motions  which  subsequently  arise,  we  have 
studied  certain  time  history  responses  of  aircraft  H  in  order  to  see 
whether  responses  become  multiperiodic  and  tend  towards  nonperiodicity 
or  chaos  (Ruelle  (1977)).  It  was  shown  previously  (Mehra  et  al .  (1977)) 
that  most  of  the  aircraft  H  time  histories  exhibit  one  of  the  following 
two  types  of  motions: 

1)  the  state  variables  behave  in  a  noticeably  periodic  manner,  with 
an  amplitude  growth/decay  time  constant  which  is  much  greater  than 
the  interval  of  interest  over  which  the  motion  is  observed; 

2)  the  motion  decays  to  a  steady-state  equilibrium  in  which  all  of 
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the  state  variables  arrive  at  a  unique  value,  without  oscilla¬ 
tions. 

However,  there  are  initial  conditions  and  control  settings  which 
produce  a  response  that  appears  to  be  erratic  and  it  is  not  possible 
to  conclusively  determine  whether  or  not  the  motion  possesses  limit  cycle 
behavior,  based  purely  on  an  inspection  of  time  histories.  Ruelle  (1977) 
states  that  when  a  motion  contains  at  least  three  basic  frequencies,  it  is 
possible  for  the  periodic  response  to  appear  random  in  nature,  due  to  cer¬ 
tain  nonlinear  perturbations. 

It  is  therefore  worthwhile  to  look  at  the  spectra  of  these  responses. 
Aircraft  H  has  been  used  in  this  study,  because  it  models  adequately  the 
kind  of  nonlinear  terms  which  generate  periodic  behavior.  The  spectra 
used  here  are  merely  the  Fourier  transforms  of  the  autocorrelation  func¬ 
tions  of  each  of  the  state  variables  (p,q,r,a,B)  computed  from  their  time 
histories.  A  routine  from  the  IMSL  library  package  was  used  to  generate 
the  spectra. 

The  control  values  for  which  spectra  of  resulting  time  histories  were 
computed  are  shown  in  Fig.  4.1.  These  figures  are  taken  from  the  report 
by  Mehra  et  al.  (1977). 

Fig.  4.2  shows  a  time  history  which  is  well-behaved.  The  small  tran¬ 
sient  behavior  at  the  beginning  decays  rapidly  to  steady  state  values  for 
each  of  the  states.  The  controls  for  this  case  are  held  at  6a  =  0°,  6e  =  2°, 
6r  =  5° ,  and  the  initial  conditions  are  a  =  8  =  p  =  q  =  r  =  0.  The  corresponding 
spectra  are  shown  in  Fig.  4.3.  These  are  generated  for  a  time  interval  from 
2  to  100  seconds  which  includes  the  initial  transient  motions.  Note  that 
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there  are  2  frequencies,  a  secondary  one  which  has  roughly  twice  the  value 
of  the  fundamental  frequency,  the  latter  being  about  0.40  cycles  per  sec¬ 
ond  (cps).  Note  also  that  the  secondary  frequency  appears  only  in  the 
spectra  of  the  longitudinal  variables,  a  and  q.  The  control  settings  for 
this  case  are  in  the  linear  region,  so  such  decoupling  is  not  a  sur¬ 
prise.  Averaging  effects  have  reduced  the  amplitude  of  the  spectra  (i.e., 
for  most  of  the  time  period  over  which  the  spectra  are  computed,  the  system 
is  at  equilibrium.  The  amplitude  thus  changes  with  the  ratio  of  transient 
time  interval  to  total  timer  interval ) . 

Fig.  4.4  shows  the  results  of  computing  the  spectra  over  the  time  period 
22  to  85  sec.  which  effectively  avoids  all  transient  behavior.  The  mag¬ 
nitudes  have  been  reduced  several  orders  of  magnitude,  and  the  secondary 
frequency  has  vanished.  Note,  however,  that  the  fundamental  frequency 
is  still  detectable,  and  it  still  has  the  same  value  (about  0.40  cps). 

Limit  cycle  motions  were  next  studied  by  repeating  the  above  run, 
except  for  6a  which  is  set  at  -18°.  Since  the  initial  conditions 
are  again  at  the  origin,  a  Hopf  Bifurcation  (i.e.,  jump)  to  a  limit 
cycle  occurs  for  this  value  of  6a,  as. seen  in  Fig.  4.5.  The  resulting 
spectra,  shown  in  Fig.  4.6,  indicate  that  the  fundamental  frequency 
is  about  1.3  cps,  with  a  secondary  frequency  at  about  2.3  cps  and 
hints  of  a  third  frequency  (see  the  p  spectrum)  at  1.05  cps.  The  major 
difference  from  the  6a =  0°  case,  however,  lies  in  the  amplitudes  of  the 
spikes,  which  are  at  least  three  orders  of  magnitude  greater.  The  spikes 
are  sharper,  indicating  a  "clean"  oscillation  (the  spectrum  of  a 
pure  sinusoid  is  an  impulse  located  at  the  frequency  of  the  sinusoid;  the 


spectrum  of  white  noise  is  a  constant  over  all  frequencies  and  is  propor¬ 
tional  to  the  power  of  the  signal) 

A  run  was  next  made  with  6a  =  -6°,  6e  =  2c',  ->  =  0%  and  initial  conditions 
p  -  -46.3,  q  =  23.4,  r  =  59.2  deg/sec,  i  =  -2.39°,  .-  =  -31.5°.  The  time  history. 
Fig.  4.7,  indicates  evidence  of  at  least  2  frequencies.  The  corresponding 
spectra.  Fig.  4.8,  show  a  dominant  frequency  at  1.2  cps,  and  three  sub¬ 
harmonics  at  about  0.15,  0.3  and  0.45  cps.  As  in  all  cases  except  Fig.  4.4, 
the  time  interval  for  the  spectrum  computation  in  2  to  100  sec.  The  some¬ 
what  irregular  spacing  of  the  high  peaks  in  Fig.  4.7  gives  rise  to  the 
cluster  of  3  secondary  peaks  in  the  spectra.  The  absence  of  rudder  in 
this  case  may  explain  the  lack  of  secondary  frequencies  in  the  spectra  for 
B  and  r,  the  two  variables  most  directly  dependent  on  <s r . 

Fig.  4.9  shows  a  more  erratic  motion,  and  its  spectrum  is  in  Fig.  4.10. 
The  initial  values  for  the  state  variables  are  again  zero,  and  6a  =  0°, 

6e  =  2°,  and  6r  is  moved  to  10°.  Again  the  roll  motion  (p)  is  dominant 
with  a  primary  frequency  of  0.375  cps.  The  pitch  rate  (q)  motion  is  less 
periodic  and  has  several  frequencies.  It  is  not  clear,  however,  whether 
any  chaotic  regimes  exist  for  aircraft  H.  McLaughlin  and  Martin  (1974, 

197b)  show  that  in  fluid  flow,  chaotic  motions  can  result  either  via  a 
phenomenon  known  as  an  inverted  Hopf  bifurcation  (i.e.  ,  existence  of  un¬ 
stable  limit  cycles  at  control  values  below  the  control  values  for  which 
a  pair  of  complex  conjugate  eigenvalues  crosses  the  imaginary  axis)  or  via 
normal  Hopf  bifurcation  exceeding  three  in  number.  (Ruelle  (1977)  has 
shown,  however,  that  chaotic  motions  or  "strange  attractors"  are  possible 
even  after  three  Hopf-type  bifurcations.)  Since  the  order  of  the  aircraft  H 
model  is  five  and  at  most  two  pairs  of  complex  conjugate  eigenvalues  can 
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cross  the  imaginary  axis,  "strange  attractors"  cannot  be  present  if 
only  normal  bifurcations  are  considered.  At  this  time,  it  does  not 
seem  that  inverted  bifurcations  are  present  in  the  aircraft  H  model, 
though  this  does  require  further  study. 
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Fig.  4.7:  A/C  H  Time  History;  <5a  =  -6°,  6e  »  2°,  6r  ■  0°; 

p  =  -46.3,  q  -  23.4,  r  -  59.2  deg/sec;  a  *  -2.398,  6  -  -31.5°. 
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Fig.  4.9:  A/C  H  Tijne  History;  6a  »=  0°,  6e  ■  2°,  6r  *  10°; 
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initial  conditions  at  origin. 
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4. 2  Command  and  Stability  Augmentation 

Progress  of  a  significant  nature  in  this  area  awaits  the  development 
of  the  spin  bifurcation  algorithm,  discussed  in  Sec.  2.1.  Using  the 
current  algorithm  several  computer  runs  are  required  in  order  to  generate 
adequately  all  of  the  branches  of  the  bifurcation  surfaces  in  the  two- 
dimensional  control  space.  When  the  third  control  is  varied,  an  entirely 
new  surface  must  be  generated. 

However,  the  results  developed  so  far  for  two-dimensional  cormiand 
augmentation  systems  using  the  aircraft  H  model  of  Mehra,  Kessel  and 
Carroll  (1977)  have  been  very  promising.  The  (two-dimensional)  bifur¬ 
cation  surfaces  are  used  to  define  relationships  between  the  two  control 
variables  which  serve  to  expand  as  much  as  possible  the  region  in  the 
equilibrium-state  space  for  which  bifurcations  are  avoided.  Work  has 
been  centered  to  date  on  the  control  pair  (5a,  5r).  Command  augmentation 
gains  relating  6r  to  5a  are  generally  called  aileron-rudder  interconnect 
(ARI)  systems.  API  gains  cause  the  rudder  to  deflect  in  conjunction  with 
aileron  movement.  The  purpose  is  to  compensate  in  some  manner  for  the 
effect  of  changing  flight  conditions  on  control  response  of  the  aircraft. 

The  standard  method  of  defining  the  ARI  gains  is  to  set  them  as  linear 
functions  of  angle-of-attack.  By  using  BACTM,  in  particular  the  bifurcation 
surface  plots,  it  is  possible  to  generate  directly  ARI,  or  any  other  type 
of  command  augmentation,  functions.  This  is  a  more  general,  or  global, 
procedure  than  those  which  reduce  to  a  gain  linear  in  a.  It  does  not 
rely  on  localized  analyses  throughout  the  control -state  space,  but  incor¬ 
porates  the  global  aircraft  behavior  information  inherent  in  the  BACTM  results. 
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Work  done  to  date  has  used  aircraft  H  bifurcation  surfaces.  For 
each  value  of  6e,  a  surface  was  generated,  and  a  linear  relation¬ 
ship  between  rudder  and  aileron  was  derived.*  The  criterion  was 
to  expand  as  much  as  possible  the  "non-catastrophic,  non-limit-cycle" 
region  in  the  control  space.  Since  each  setting  of  6e  corresponds  to  an 
equilibrium  value  of  angl e-of-attack,  a,  the  BACTM  ARI  gains  can  similarly 
be  plotted  versus  a.  Fig.  4.11  shows  this  plot  compared  to  the  linear  ARI 
gains  selected  by  Gilbert,  Nguyen  and  VanGrunst  (1976).  The  main  point 
to  be  made  here  is  that  the  general  sense  of  the  two  plots  is  sinilar. 

(In  Fig.  4.11,  the  values  at  the  break  points  refer  to  the  elevator  deflec¬ 
tion,  in  degrees).  The  gain  values  are  of  comparable  magnitude,  and  a 
"mean  slope"  for  the  BACTM  points  would  not  be  very  different.  It  is 
felt  that  the  BACTF1  method  would  result  in  better  overall  performance 
because  of  the  global  stability  information  inherent  in  it.  Work  will 
be  continued  in  this  area  with  the  F-4  model. 


*BACTM  does  not  require  that  this  relationship  be  linear.  Later  results 
may  show  better  response  for  gains  nonlinear  in  (6a,6r). 


CHAPTER  V 


Conclusions  and  Recommendations  for  Further  Research 

5. 1  Conclusions 

Based  on  analysis  using  BACTM  of  the  aircraft  F  model  defined  in 
this  report,  we  can  conclude  that 

(i)  BACTM  has  been  succesfully  expanded  and  modified  so  that  it 
now  has  the  capability  to  perform  analysis  of  aircraft  in  the 
highly  complex  post-stall  and  spin  regimes.  This  was  achieved 
by  introducing  more  powerful  continuation  methods  for  solving 
equilibrium  and  bifurcation  surfaces,  and  by  utilizing  an 
aircraft  model  with  sufficient  aerodynamic  data  to  simulate 
motions  over  extreme  ranges  of  angle-of-attack  and  sideslip. 

(ii)  The  aircraft  F  model  was  very  useful  for  its  role  in  the  de¬ 
velopment  of  the  BACTM  spin  analysis  program.  Also,  it  was 
useful  for  studying  developed  spin  motion.  However,  the  aero¬ 
dynamic  data  as  given  cannot  with  uniform  accuracy  deal  with 
the  wide  variation  in  flight  conditions  which  results  from 
maneuvers  proceeding  from  trim  to  spin  conditions.  It  is 
necessary  to  begin  using  a  more  complete  model,  with  the 
aerodynamics  divided  into  static,  forced  oscillation  and 
rotary  balance  data.  In  this  way,  combinations  of  the  three 
sets  may  be  varied  from  one  flight  regime  to  another,  to  more 
accurately  simulate  actual  flight  test  results. 

(iii)  Spin  study  via  BACTM  can  be  made  easier  by  assuming,  on  a 
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first-trial  basis,  that  gravity  effects  are  negligible. 

This  assumption  must  be  applied  with  care,  however,  because 
gravity  plays  an  important  role  in  post-stall  dynamics  as 
well  as  most  spin  motions.  Nonetheless,  it  appears  to  be 
reasonable  to  apply  this  assumption  to  aircraft  F  studied  here. 

(iv)  The  transition  between  the  non-spin  and  spin  stable  equilibria 
for  aircraft  F  is  difficult  to  achieve  due  to  the  strongly  attrac¬ 
ting  nature  of  the  intermediate,  high-a  equilibrium  region.  The 
BACTM  analysis  shows  a  high  degree  of  nonlinear,  oscillatory,  limit 
cycle  behavior  associated  with  a  large  domain  of  attraction  and  a 
large  region  of  structural  stability  for  the  limit  cycle  family. 

(v)  Aircraft  F  in  flat  developed  spin  follows  a  tight  vertical 
helical  path,  which  is  characterized  by  constant  speed, 
sink  rate,  and  a  high,  autorotational ,  steady  yaw  rate. 

In  this  condition,  the  elevator  is  ineffective  as  a 
control,  and  recovery  is  possible  only  by  using  the  rudder 
and  aileron. 

(vi)  Entry  into  spin  for  aircraft  F  is  strongly  affected  by  assump¬ 
tions  about  velocity.  A  spin  of  smaller  amplitude  oscillations 
and  higher  angle-of-attack  (i.e.,  a  more  flat,  developed  spin) 
results  when  velocity  is  fixed,  as  opposed  to  when  it  is  allowed 
to  vary.  The  fixed  velocity  case,  for  small  a  and  b 

angles,  corresponds  to  a  nonzero  thrust  situation,  with  the 
thrust  neutralizing  drag  effects.  When  velocity  is  allowed 
to  vary  in  this  report,  zero  thrust  is  assumed. 


251 


(vii)  The  stall  departure  region  is  preceded  by  a  wing  rock 
type  of  limit  cycle  motion  which  is  typical  of  modern 
high  performance  aircraft. 


5.2  Recommendations 

Based  on  our  experience  during  this  reporting  period,  it  is  suggested 
that  the  following  areas  be  investigated  in  the  future: 

(i)  The  F-4  data  should  be  analyzed  using  BACTM  since  it  contains 
more  realistic  aerodynamics,  is  well  documented  and  well 
supported  by  flight  tests. 

(ii)  Computational  development  of  BACTM  should  be  continued,  as 
more  accurate  and  efficient  algorithms  are  needed,  particu¬ 
larly  for  the  generation  of  a  full  set  of  bifurcation 
surfaces. 

(iii)  More  analysis  should  be  made  of  the  assumptions  regarding 
equilibrium  motion  in  spin,  particularly  the  assumption  that 
gravity  can  be  neglected  in  certain  cases. 

(iv)  More  time  history  runs  and  analyses  should  be  done  with  initial 
conditions  in  selected  regions  of  the  state-control  space; 

the  purpose  being  to  more  clearly  define  persistence  of  limit 
cycle  behavior  and  to  establish  boundaries  for  the  domains  of 
attraction.  Analytical  and  computational  procedures  need  to 
be  developed  for  generating  Hopf  bifurcation  surfaces. 
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(v)  Using  aircraft  F-4  data  to  generate  bifurcation  surfaces, 
perform  preliminary  synthesis  of  a  command  augmentation 
system  using  BACTM.  This  system  should  then  be  compared 
to  other  systems  in  the  literature. 

(vi)  Using  the  F-4  model,  study  the  role  of  thrust  in  post-stall, 
departure,  spin  entry  and  developed  spin  flight  conditions. 

(vii)  Determine  the  parameter  values  under  which  structurally  stable 
limit  cycles  such  as  high-a  oscillatory  spins  exist  for  a 
given  aircraft  model  and  design  dynamic  control  strategies 
for  recovery  from  such  limit  cycles. 


! 


»o 


APPENDIX  A 


Notation 


b 


M 

C  =  — 
l  qSb 

M 

C  =  -X 
m  qSc 

M 

C  =  -i- 

qSb 

(F-mg)-x 
C  =  ~ 

x  qS 

(F-mg)*y 
r  =  ~ 

y  qs 

(F-mg)*z 
Cz  =  qs 

E 

f 

F 

F 

9 

9 

G 


wing  span 

mean  aerodynamic  chord 
rolling  moment  coefficient 

pitching  moment  coefficient 

yawing  moment  coefficient 

longitudinal  force  coefficient 

side  force  coefficient 

normal  force  coefficient 

vehicle  total  kinetic  energy 
force-moment  terms  in  the  aircraft 

Jacobian  matrix  of  partial  derivatives, 

j 

aerodynamic  force 

(constant)  acceleration  due  to  gravity, 

9.8067  m/sec2  (32.174  ft/sec2) 

algebraic  system  of  terms  for  generating  bifur¬ 
cation  surfaces  (Chapter  II) 

augmented  Jacobian  matrix  of  partial  derivatives, 

r*,l 

JjjyjJ  »  for  bifurcation  surfaces 
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W  !z  ’ 1  xz 


Mx’My’Mz 


m 

M 

N(-) 


P,q»r 


q 

R 

R(-) 


RjTjZ 


I 


s 

t 

t* 


u,v,w 


V 

Vl 

V 

w 

x 


altitude  above  earth's  surface 

body  axis  moments  and  product  of  inertia,  taken 
about  the  center  of  mass 

moment  of  inertia  tensor  (Eq.  3.2.38) 

rolling,  pitching,  yawing  moments  acting  about  body 
axes 

aircraft  mass 

special  vector.  Sec.  2.1.2 
null  space 

angular  rates  about  body  axes  (roll,  pitch,  yaw, 
respectively) 

dynamic  pressure,  %pVz 

radius  of  helical  path  of  airplane 

range  space 

n-dimensional  space  of  real  numbers 

unit  vectors  in  cylindrical  coordinates;  z1  is 
vertical,  directed  toward  center  of  earth" 

wing  area 

time 

time  at  which  equilibrium  solution  is  made 
body  axis  components  of  V 
airspeed,  =  |V| 

horizontal  component  of  velocity 

aircraft  center  of  mass  velocity,  inertial  with 
respect  to  local  horizontal 

aircraft  weight 

vector  of  state  variables;  e.g. ,  for  na  5  equili¬ 
brium  system,  x =  (p,q,r,a,6) 
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a 


x,y,z 


a 


8 

5 
A 

6 

6a, 6e,6r 


P 


ip,d,<p 


0) 

n 


a  e  [-1,1] 
a  e  (-1,1) 


augmented  vector  of  dependent  variables,  for 
bifurcation  surfaces;  e.g. ,  for  n  =  5  equilibrium 
system,  y  =  (x,6. ),  where  6.  e  (6 a,6e,6r) 

”  J  J 

aircraft  body  axis  unit  vectors  (x  positive  through 
nose,  y  positive  through  right  wing,  z  positive 
down ) 

angle  of  attack,  or  incidence  angle  (Chapter  III); 
also  continuation  variable  (Chapter  II) 

angle  of  sidesl  ip 

control  parameter;  either  6a,  6e,  or  6r  (Chapter  II) 
determinant  of  F,  the  Jacobian  matrix 
control  vector,  (5a,6e,6r) 

aileron,  elevator,  and  rudder  control  deflections 
(positive  6e  is  trailing  edge  down,  positive  6a  is 
right  trailing  edge  down,  positive  6r  is  trailing 
edge  left) 

atmospheric  density 

Euler  angles  defining  orientation  of  body  axes 
in  the  inertial  reference  axes  (yaw,  pitch,  roll, 
in  that  sequence) 

angular  rate  about  center  of  mass,  /pz  +  q*  +  ri 

polar  angle  in  cylindrical  coordinate  system  defining 
aircraft  position 

a  such  that  -1 <  a <  1 

r  such  that  -1  <  a  <  1 


the  combinations  [•)  and  (•]  are  similar;  i.e.,  ae[-l,l)  means 
-1  <  a  <  1 


a  e  A  a  is  an  element  of  the  set  A 

a  i  A  a  is  not  an  element  of  the  set  A 

[a.  .]  a  matrix  array  whose  elements  are  designated  by 

a.,,  the  element  at  row  i,  column  j 

*  sJ 
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(a,) 


V 


I 


( ) 


T 


(  )* 


det(* ) 


C) 

& 


a  vector  array  whose  elements  are  designated  by 
a..,  the  i  th  location  element 

indicates  vector  v  is  in  inertial  coordinates 

matrix  transpose 

i)  complex  conjugate,  as  in  Eq.  (2.1.41) 

ii)  equilibrium  solution,  as  in  Sec.  2. 1.4.1 

the  Euclidean  norm  of  the  vector  x,  i.e., 

II  *  II  =J I  if  ?eRn 

Vi  =  l 

the  determinant  of  the  argument  (which  must  be  a 
square  matrix) 

d(  )/dt;  also  d(  )/ds  in  Sec.  2.1.1 
equal  by  definition 


Stability  Derivatives 


A  3Ci 

C1  -  ,  for  i  =  £,m,n,x,y,z 

and  £  =  <5a,6e,<5r 

A  3Ci 

C  Q  — K —  ,  for  n  =  P.r,§ 

’n  3(^n) 


ap 

Ci  -  — t —  ,  for  s  =  q 
^  3(^0 


In  addition,  the  coefficients  C.  and  the  derivatives  are 

i  ' 

functions  of  a  and  6.  and  are  usually  determined  via  tabular  or  graphical 
look-up. 
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